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1. INTRODUCTION 


The General Interpolants Method (GIM) computer code was developed 
to analyze complex flow fields which defy solution by approximate methods. 
The code uses numerical techniques to solve the full three-dimensional time- 
averaged elliptic Navier-Stokes equations in arbitrary geometric domains. 
The equations are cast in Conservation Law Form and written in an orthog- 
onal Cartesian coordinate system. A generalized three-dimensional geom- 
etry package is used to model the flow domain, generate the numerical grid 
of discrete points and to compute the local transformation metrics. Com- 
putation is done in physical space by explicit finite -difference operators. 

The approach taken in developing GIM is akin to the Method of Weighted 
Residuals which can produce discrete analogs classically termed finite - 
element methods. The GIM approach essentially combines the finite ele- 
ment geometric point of departure with finite difference explicit computa- 
tion analogs. This provides a capability which takes advantage of the 
geometric flexibility on an element description and the superior computa- 
tion speed of difference representations. 

The numerical analogs of the differential equations are derived by 
representing each flow variable with general interpolation functions. The 
point of departure then requires that a weighted integral of interpolants be 
zero over the flow domain. By choosing the weight functions to be the 
interpolants themselves, the GIM formulation produces identically the clas- 
sical implicit finite-element discrete equations. We do not use these forms 
in the GIM/STAR code due to their fully implicit nature. Rather, we choose 
the weight functions to be orthogonal to the interpolant functions. This 
choice reproduces explicit finite -difference type discrete analogs. By ap- 
propriate choice of constants in the orthogonal weight functions, the GIM 
becomes analogous to known finite difference schemes such as centered 
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differences, windward schemes, and the two-step MacCormack technique. 

The GIM analogs, however, are automatically produced for arbitrary geom- 
etric flow domains. For pure Cartesian domains, the GIM method exactly 
reproduces these known schemes, but is a more general point of departure 
and provides greater flexibility in choice of difference scheme. 

A motivation for developing this code on these principles was to 
provide an analytical tool which is more user oriented than the basic re- 
search tools which exist. A fully production-line code to solve the com- 
plex Navier -Stokes equations does not exist today. In developing the GIM 
code, we have attempted to bridge the gap somewhat between the pure re- 
search codes and the ultimate production tool. The code was originally 
developed for a CDC 7600 computer system. It has subsequently been con- 
verted to execute on Univac 1108 and 1110 systems. The most recent ac- 
tivity, which is the subject of this contract, has been the conversion of the 
GIM code to the STAR- 100 computer at NASA-Langley Research Center. 

Development of the GIM/STAR code is motivated by a number of factors: 
(1) NASA -Langley currently has a need for an analysis tool to aid in its re- 
search programs in hypersonic aerodynamics and propulsion; (2) The STAR 
machine has the potential of increased speed over a 7600 for flow field com- 
putation; (3) the virtual memory concept of the STAR machine is well suited 
for codes requiring large data bases; and (4) the GIM formulation and con- 
cepts are readily vector izable for a machine such as the STAR. A feasibility 
study of the compatibility of the GIM concept with the STAR machine resulted 
in the conclusion that definite and significant advantages could be achieved. 

The objectives of this contract were threefold; (1) convert the GIM 
code into STAR vector FORTRAN and implement it on the computer system 
at NASA -Langley; (2) upgrade the code to include a two- gas capability and 
a sharp corner boundary treatment; and (3) provide a compatible graphics 
package for operation on the CDC 6000 series machines. These objectives 
have been met with the code being successfully implemented on the STAR 
system. Two primary example problems have been computed; (1) a 


1-2 



two-dimensional simulation of a Scramjet exhaust interacting with a super- 
sonic free stream flow, and (2) a three-dimensional Scramjet simulation. 
Comparison of the STAR CP run times for these and other problems with 
a CDC 7600 have given improvement factors of up to 5 for the STAR code. 

The resulting calculations are identical to the 7600 solutions. The wall- 
clock time required to complete a GIM/STAR run is also significantly lower 
than a typical 7600 job for cases which do not require large page faulting on. 
STAR. Large page faults occur when a problem size is too large to fit "into 
core" and must be stored on external media. For larger problems which re- 
quire page faults, the throughput time is comparable to the GIM/7600 code. 

This ratio could be significantly improved by adding more high speed memory 
to the STAR system. The additional technical capability tasks and the graphics 
module addition have also been implemented and verified. 

This report documents the results of the contract and also serves as 
a User's manual for the GIM/STAR code. Decks for the program modules 
reside on permanent files in the STAR computer system. FORTRAN listings 
can be obtained from these files. Section 2 of this report summarizes the 
contract study results and conclusions. Section 3 presents a general dis- 
cussion of the GIM formulation, the structure of the code for STAR and ex- 
ample runstreams for executing the code. The input guides constitute Sections 
4, 5 and 6 for the Geometry, Integration and Graphics Modules respectively. 
Hints on using the code and example problem setups are given in Sections 7 
and 8. Section 9 contains a list of referenced documents which provide more 
details of the GIM theory and methods development. 

The version of the GIM code documented in this report is designated 
SE-1 (STAR -Elliptic No. 1). The User's Manual is valid only for this version 
of the GIM code. 
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2. SUMMARY OF CONTRACT STUDY 


2.1 OBJECTIVES 

The objective of this contract is threefold: (1) code the GIM formulation 
to allow longer vector lengths and to take full advantage of the STAR vector- 
ized multiplications; (Z) update the GIM code to allow the treatment of internal 
sharp corners and also to include a two-gas capability (each having different 
transport properties) for computing exhaust/freestream interaction; (3) de- 
liver the modified code and a compatible graphics package (to operate on CDC 
6000) with appropriate documentation and demonstrate its capabilities through 
the use of a sample case prepared by NASA -Langley. 

The following tasks were performed to accomplish these objectives: 
Task 1 


The standard Lockheed-Huntsville elliptic GIM code was implemented 
on the STAR system. This capability includes arbitrary two- and three- 
dimensional geometries, inviscid Euler and full Navier-Stokes equations, 
laminar viscosity model, ideal gas law for two streams with different specific 
heat ratios and explicit finite -difference analogs. Four modules were imple- 
mented; (1) module 1 — geometry and mesh generation; (2) module 2 — matrix 
assembly module; (3) module 3 — unsteady integration module and (4) Module 
4 — graphics package on CDC 6000. 

Task 2 


The current standard GIM code was upgraded to handle gas stream 
interactions with different specific heat ratios, y. An additional equation 
was added to compute the relative composition of the two gases at local 
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grid stations. The ideal gas law is used to relate the local gamma in the 
shear layer to the two "pure" gas gammas, and compute pressure from die 
energy, density and velocity fields. The second modification that was made 
is an option on handling flows over sharp corners through the addition of a 
Prandtl-Meyer like expansion equation. 

Task 3 

A graphics package compatible with GIM outputs was implemented on 
the Langley CDC 6000. 

Task 4 


The GIM code as implemented on the STAR computer and associated 
graj^ics package as implemented on CDC 6000 is documented in a User's 
Manual and Input Guide. 

Task 5 

Test cases (two) which were computed by the GIM code on the STAR 
computer system are two- and three-dimensional Scramjet exhaust plume 
problems already analyzed with GIM on a CDC 7600 as reported in Ref. 1. 
The sample cases were compared with ones previously run on a CDC 7600 
system for a comparison of run time and computational results. 

2.2 RESULTS 

The GIM code was implemented on the STAR- 100 system at Langley. 
The following are the major modifications which were made to the code; 

• Combined GEOM and MATRIX Modules to Eliminate Large 
Data File Transfers. 

• Recoded Portions of GEOM Module to Use STAR Vector 
FORTRAN and Interleaved Large Pages. 

• Added Algorithms to Diagonalize Coefficient Matrices into 
One -Dimensional Arrays. 
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• Eliminated Implicit Solution Algorithms {Finite Element). 

• Recoded INTEG Module in STAR FORTRAN. Eliminated Large 
Data Files — Used Large Core. 

• Added Algorithms to Use Operations with Vector Lengths Equal 
to Number of Nodes. 

• Restructured Large Arrays in Common Blocks to Interleave 
Data Base Layout. 

Most of these modifications were programming changes. There was, how- 
ever, one item which was completely reformulated. This item is the matrix 
storage and multiplication operations. Scalar operations on the 7600 code 
were replaced with vector multiplies and adds to take full advantage of the 
STAR features. This basic operation involves the following calculations: 
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Most of the CP time is spent in forming the vector products and sums, i.e.. 



The GIM/STAR code was vectorized by arranging the B, C,D matrices as 
shown in Fig. 2-1. The order of operation proceeds down each of the 
diagonals with elements a,b, c,d, respectively. The elements designated 
with X are off-diagonal terms due to boundary conditions. These are 
treated by scalar operations. The example shown in this figure is for a 
two-dimensional, 20-node sample case. Three-dimensional cases can 
have up to eight diagonals. 
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Fig. 2-1 - Matrix Diagonalization Technique Used in Vectorizing 
the Basic GIM Code Calculation 
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The vectorized GIM/STAR was run using four test problems. The 
first two of these are "simple" cases used mostly for debugging. They 
consist of two-dimensional supersonic source flow in a 15 deg. duct and a 
converging-diverging nozzle flow. The two primary test cases consist of 
the problems already solved with GIM and reported in Ref. 1. These are 
termed the two-dimensional Scramjet and the three-dimensional Scramjet 
cases* respectively. The GIM/STAR calculations reproduced the results 
previously obtained on the CDC 7600. 

A study was made of CP time improvement of the GIM/STAR code over 
the 7600 version. The results of this study are shown in Fig. 2-2 for the four 
test cases just described. This data is shown in Fig. 2>3 plotted as CP times 
versus number of nodes in the problem. As seen from these data, a CP ratio 
of 5.5 was obtained for test case 1 with ratios generally from 3 to 5 for the 
other problems. This is a significant improvement in CP time especially 
since GIM/STAR is coded entirely in vector FORTRAN. 




CP Seconds for 
100 Time Steps 


Case 

Nodes 

7600 

STAR 

Ratio 

2-D Source Flow 

121 

4.2 

2.0 

2.1 


442 

15.4 

4.0 

3.9 


961 

33.6 

6.8 

4.9 


1681 

58.8 

10.6 

5.5 

2-D Nozzle 

1932 

70.8 

13.6 

5.2 

2-D Scramjet 

940 

29.3 

9.6 

m 

3-D Scramjet 

7904 

480.0 

124.0 

3.9 


Fig. 2-2 - Summary of CP Time Comparison of GIM/STAR Code 
and CDC 7600 
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CP Seconds (x 10 ) per Node-Step 



Fig. 2-3 - CP Time vs Problem Size for STAR and 7600 Computer 


The other significant study of machine efficiency involves the "large 
page fault" feature of the STAR virtual memory system. For problems 
which are small enough to fit into "core," i.e., no large page faulting, the 
throughput time of a GIM/STAR job is approximately two times faster than 
a typical CDC 7600 job. For large problems which require some large page 
faults, the throughput time is about the same as a typical CDC 7600 job. The 
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cost of a GIM/STAR job is then strongly influenced by the size of problem, 
i.e., number of nodes. 


2.3 CONCLUSIONS AND RECOMMENDATIONS 

All tasks of the contract have been completed. This report constitutes 
the final report and User's Manual as required by the contract. The follow- 
ing is a summary of the conclusions reached and recommendations for further 
utilization of the GIM code on the STAR computer at NASA -Langley. 

• The GIM formulation proved to be ideally suited to the structure 
of the STAR -100 system. 

• Vectorization was easily implemented using vector lengths equal 
to the number of nodes. 

• STAR FORTRAN coding resulted in CP time ratios of approxi- 
mately 5 over a CDC 7600 version of the code. 

• Problem sizes which do not require large page faults run about 
two times faster than a typical 7600 GIM run. 

• Large problem sizes do cause page faulting of the GIM/STAR 
code and hence total run times comparable to the 7600 code. 

• The GIM/STAR code has successfully computed Scramjet nozzle- 
exhaust flows which have not been previously attempted. 

• For flows which are entirely supersonic, a hyperbolic GIM algo- 
rithm should be added to the STAR code to reduce the cost of a 
run. 

• For viscous flow in which downstream influences are negligible, 
a GIM parabolic solution algorithm should be added to further 
reduce the run time. 

• A "partially-elliptic" algorithm can also be added to solve vis- 
cous flows with some downstream feed-back effects. 

• The special STAR language SLl should be investigated for potential 
use in GIM to further reduce the run time. 

• Turbulence models, eddy viscosity and TKE, should be incorporated 
in GIM/STAR to allow turbulent viscous flows to be computed. 

• An equilibrium chemistry model should be added in addition to the 
ideal gas model now used. 

• The GIM/STAR code should now be run for additional problems such 
as Scramjet inlet flows, flows in the combustor area and some sim- 
plified external aerodynamic flows. 
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3. GIM CODE STRUCTURE 


3. 1 SUMMARY OF GIM THEORY 

A formal development of the General Interpolants Method (GIM) is 
given in Refs. 1 and 2. The following discussion summarizes the theory, 
equation sets and numerical approach used in the GIM code for use on the 
STAR- 100 computer. 

GIM is a new methodology for constructing numerical analogs of the 
partial differential equations of continuum mechanics. A general formula- 
tion is provided which permits classical finite element methods and many 
of the finite difference methods to be derived directly. The GIM approach 
is new in the sense that it combines the best features of finite element and 
finite difference methods. The technique allows complex geometries to be 
handled in the finite element manner and operates on the integral form of 
the conservation laws. Solutions can be generated implicitly with the finite 
element analogs or by explicit finite difference analogs, which do not require 
a reduction of large systems of linear algebraic equations (no matrix inverse). 
The STAR version of GIM contains only the explicit finite difference analogs. 

As is the case with all attempts to solve partial differential equations 
by numerical approximations, the domain of interest is first discretized by 
appropriate subdivision into an assemblage of interconnected finite elements. 

A mesh generation is used in the GIM approach which incorporates general 
curvilinear coordinates, stretching transformations and bivariate blending 
to produce an automated mesh/element generation. Shape functions based 
on a set of generalized interpolants are then chosen to describe the behavior 
over each element. We then proceed, as in the Method of Weighted Residuals 
(Ref. 3), by multiplying the discretized equations by a set of weight functions 


3-1 



and integrating over the volume of the element. By choosing the weight 
functions equal to the shape functions, we reproduce via Galerkin (Ref. 3), 
the classical finite element nodal analogs. It is at this point that we intro- 
duce one of the important concepts of GIM: orthogonal weight/ shape func- 
tions. By appropriately choosing the weight functions to be orthogonal to 
the shape functions, we can obtain explicit nodal analogs. Further, by a 
choice of arbitrary constants in the orthogonal weight functions , we can 
produce known finite difference nodal analogs, such as centered difference, 
upwind/downwind differences and the two-step MacCormack algorithm 
(Ref. 4). Asa result of this spatial discretization, we have reduced the 
partial differential equations to ordinary differential equations with "time" 
as the independent variable. Any forward marching algorithm such as Euler, 
Runge-Kutta or predictor-corrector can be used to advance the solution pro- 
files in time. 

The GIM formulation is not a Finite Element method in the classical 
sense. Rather, finite difference methods are used exclusively but the equa- 
tions are written in general orthogonal curvilinear coordinates. Transform- 
ations are used to transform the physical planes into regions of unit cubes. 

The mesh is generated on this unit cube and the local metric coefficients 
generated. Each region of the flow domain is likewise transformed and then 
blended to form the full flow domain. In order to treat "completely-artibrary" 
geometric domains, different transformations may be employed in different 
regions. For this reason, we then transform the blended domain back to phys- 
ical space. This allows the same set of equations to be solved in each region, 
with the local Jacobian of the transformation being the coefficients. 

Equation Sets; The partial differential equations solved by the GIM code 
are the Navier-Stokes written in three-dimensional Conservation Law Form for 
a Cartesian coordinate system. The full set of differential equations solved 
by the GIM code is shown in Fig. 3-1. These will be referred to as the Navier- 
Stokes equations, although they include continuity and energy conservation as 
well as momentum. For inviscid flow calculations without shock waves, the 
j terms can be set to zero to produce the Euler set. Fig. 3-1 lists the viscous 
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terms with coefficients* [it X* k for the first and second coefficients of vis- 
cosity and the thermal conductivity, respectively. These are specified by 
the user as "laminar" constants. For reasons of numerical stability and 
capture of strong shock waves, additional components of the diffusion coef- 
ficients are added automatically by the GIM code. Figure 3-2 is a list of 
currently used Numerical Diffusion Cancellation (NDC) coefficients. These 
are added to the real diffusion coefficients. The purpose of these NDC coef- 
ficients is to cancel low order truncation error terms which arise in the 
numerics. Reference 5 presents the basic principles of the NDC technique 
in more detail. 

The differential equations are solved in strong conservation or diver- 
gence law form. The solution is started at some time t where the entire 
flowfield mesh is specified. The unsteady or relaxation of the equations is 
then done using the user-specified nodal analogs. At this time, the pressure 
is computed from the ideal gas law for a single or two component gas. 

Geometry and Coordinate Tregttmentt The underlying concept of the 
theory of finite elements relates to the classical notion of "piecewise approxi- 
mation." A real continuum, R, can be considered as an infinite number of 
interconnecting points. A piecewise approximation to the continuous region, 

R, consists of a finite number of interconnecting points which, in some man- 
ner, resembles R itself. We are also concerned with a continuous field, f, 
whose domain of interest is the region R; for example, the region R can be 
a body over which a flow field is to be determined. The function f(X) pre- 
scribes a unique value of f at an infinite number of points X in R. The 
problem of piecewise approximation is to construct a discrete model of the 
region R and the function f such that the field is represented approximately 
at a finite number of points in R. 

An important concept of the theory of finite elements is the "disconnect- 
ness" property. This means that a domain can be divided into a- finite number 
of pieces called elements and the approximation to the functional distribution 


3-4 




CoeffLcients for Stress Terms 
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= laminar, constant viscosity 
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k^ = thermal conductivity 
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NDC Values 


Fig. 3>2 - NDC Coefficients Used in the GXM/STAR Code 
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over each element can be studied independently. Thus the approximating 
functions for each element completely define the behavior of the function 
profile within the element without consideration of Its ultimate location in 
the full model. After each element is defined, the complete discrete model 
of the body is obtained by "assembling the system." The assembly is per- 
formed by means of mathematical blending of each subdomain while main- 
taining continuity at the junctions. 

The General Inter polants Method uses these concepts, borrowed from 
finite element theory, to obtain discrete finite difference models of the Navicr- 
Stokes equations in arbitrary geometric domains. The developnient is done 
in local curvilinear intr insic coordinates based on; 

1. Analytical regions such as rectangles, spheres, cylinders, hexa- 
hedrals, etc., which have intrinsic or natural coordinates. 

2. Complex regions that can be subdivided into a number of smaller 
regions which can be described by analytic functions. The de- 
generate case is to subdivide small enough to use very small 
straight-line segments. 

3. Intrinsic curvilinear coordinate systems that result in constant 
coordinate lines throughout a simply connected, bounded domain 
R in Euclidean space. 

4. The intersection of the lines of constant coordinates that produce 
nodal points evenly spaced in the domain R. 

5. Intrinsic curvilinear coordinate systems that can be produced by 
a univalent mapping of a unit cube onto the simply-connected 
bounded domain R. 

Thus, if a transformation can be found which will map a unit cube uni- 
valently onto a general analytical domain, then any complex region can be 
piecewise transformed and blended using general inter polants. 

Consider the general hexahedral configuration shown in Fig. 3-3. The 
local Intrinsic coordinates are t]j^, rj^, with origin at point Pj^. The shape 
of the geometry is defined by 

1. Eight corner points, through Pg as shown {Fig. 3-3a) 
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2. Twelve edge functions, E (Fig. 3- 3b) 

3. Six surface functions, S (Fig. 3-3c) 

This shape is then fully described if the edges and surfaces can be specified 
as continuous analytic vector functions 

_S.(x,y,z). Ej^(x,y,z). 

Based on the work of Gordon and Hall (Ref. 6) we have developed a general 
relationship between physical Cartesian space and local curvilinear intrinsic 
coordinates. This relation is given by the following general trilinear inter- 
polant; 

S 5 + rj 1 + (l-n^) ^2 + ^2 S 4 

+ ( 1 -H 3 ) S J + 113 

- (l-r?i) (l-rj2) E5 - (l-nj) ^2 Eg - (l-n^) 

" *^1 ^2 ^7 ” ^4 ” *^3 ^12 

- ^2 " ^^3 ^10 ’ (1-H3) E^ 

(3.1) 

- (I-T 72 ) ^3 Eg - ^ 2 ^ (l”n3) E^ " ^2 ^11 

+ (1-rjj) (l-r)2) (1-^3) Pj + (1-Hj) (1-^2) ^5 

+ (1-rij) 1I2 ^4 ^ *^2 ^^3 ^8 

+ rjj (l-r/2) (1-^3) ^2 + r»i (1-H2) ^3 
+ rij r)2 (I-H3) P3 + H2 ^3 P7 
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X (rjj, rj 2 * H 3 ) 


r / 

X (nj» M2* M3) 

y (Mj» M2» M3) 

^ (Mj» M2* M3) 


and S^, are the vector functions defining the surfaces and edges, respec- 
tively, and ,P. arc the (x, y, z) cot>rdinates of the corner points. Edge and 
surface functions are currently included in the GIM code as summarized in 
Fig. 3-4. 

• HEXAHEDRAL BUILDING BLOCKS 

• EDGES (Combinations of up to Five Types) 

Linear Segment 
Circular Arc 

Conic (Elliptical, Parabolic, Hyperbolic) Arc 
Helical Arc 
Sinusoidal Segment 

• SURFACES (Bounded by Above Edges) 

Flat Plate 
Cylindrical Surface 
Edge of Revolution 

• REGION OF INTEREST DESCRIBED BY UNLIMITED 

NUMBER OF HEXAHEDRAL BUILDING BLOCKS 


Fig. 3-4 - GIM Three-Dimensional Geometry Package 
Edge and Surface Function Types 

With tills transformation, any point in local coordinates r]^, r]^* M3 can 
be related to global Cartesian coordinates x, y, z. Likewise any derivatives 
of functions in local coordinates can be related to that derivative in physical 
space, i.e., 
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where J, the Jacobian matrix, is easily computed from Eq. (3.1); 
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Aj = det|j| 


(3.3) 


Discretization of Equations; Equations (3.2) and (3.3) then relate any 
flow gradient in t) ri space back to physical space. We then need a 
method of approximating the flowfield function f in local coordinates. This 
is done with a set of general interpolation functions I(t) j, ^ 3 ) such that 

8 

lj(n^.T) 2 »rj 3 ) f. (3.4) 

j=l ^ 


where f^ are the flow variables at the corner points of the element. The simplest 
form for the are the trivariable Lagrange inter polants. The theory itself does 
not restrict the L to be linear, but in the elliptic GIM code we currently have 
only the trilinear interpolants. Any flow gradient can then be simply computed as 

8 

9f ^ 

9n- 

j=i 


ai. 

f- 

J 


(3.5) 


The elliptic GIM code uses these concepts to produce a discrete analog 
of the Navier-Stokes equations for a single analytical region. Proceeding as 
in Ref. 1, the point of departure is the requirement that 


SR 
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W dV = 0 
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where W is an arbitrary weight function and 4^ is the differential equation 


;)U i)E ;)F (')G 

at Ox 0 ’^ Oz ■ 


(3.7) 


The general interpolants 1 are used to approximate each of the flow variables 
U,E,F,G and arc substituted into Eq. (3.6) to obtain a discrete analog; 


fA..] U. + [b..] E. + [C..] F. + [d..] G. = 0 


(3.8) 


where U., E., F., G. are the flow variables at the node points of an element 
J J J J ^ . 

and the coefficient matrices, geometrically dependent, are given by the follow- 
ing integrals. 


.1 *1 1 




(3.9) 
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The procedure is then to obtain Individual element matrices by numer- 
ical Gaussian quadrature evaluation of Eq. (3.9). We thus have a Cartesian- 
like finite difference equation (Eq. (3.8)) with all the geometric transformation 
data as coefficients of the finite difference operator. The A,B,C,D matrices 
may be viewed as ’'derivative takers" or finite difference operators which con- 
tain both the difference scheme and the transformation metrics. 

To model a complex three-dimensional donna in, it is first subdivided 
into regions and Eq. (3.8 ) obtained for each region. After all regions are 
so processed, the final order of business is to "assemble the system" to ob- 
tain a consistent sot of coefficient matrices for the full complex region. This 
is accomplished by the so-called quasi-var iational procedure. 


Since the choice of the weight function W is arbitrary, the integral, (Eq. 
(3 6)) must bo stationary for variations in W. The so-called quasi-variational 
approach determines the conditions under which the integral is stationary. 

That is to say that if W' = W + 6W then 


Now let 


.so tliat 


N ^ 

6 j 6W 4^ dx dy dz = 0 
e=l R 


W - Wj (»)j.))2.r)3) U. 


or 


6W ^ W. 6U. 

J J 




^ r 

E /'"i 


4’ dx dy dz 6U. 


= 0 


e= 1 


R 


(3.10) 


For tlie to be stationary, the sum of the coefficients of 6U^ must be 
zero. This provides a rationale for the classical assembly procedure for 
MWR derived finite element schemes. 
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In mathematics, this can be written as 

[^mn] = 2 Vi *„j (3. 1 1) 

e 

where 

1 if node L of element e coincides 
with node m of the region R 

0 otherwise 

and similarly for the B,C,D matrices. The rationale for doing this is di- 
rectly derivable from the quasi -variational procedure. 


6 . 
mi 


Upon assembly then the finite difference equations for the full domain 
have the form 


[a Iu+Tb 1e+[c 1f+ Td I 

1^ mnj n L ninj n [ innj n L mnj 


n 


= 0 


(3.12) 


At tills point, the weight functions W are arbitrary. Thus, in general, 
Eq. (3.12) is an implicit difference operator. If we choose the weight func- 
tions W to be the general interpolants I, i.e., W = I, then equation (3.12) is 
a classical Galerkin finite element model. We do not , however, use this 
description in the GIM code. Rather, we choose the weight functions W to 
be orthogonal to the general interpolants 


(W*Aj)i.I 

Now, inspection of Eq. (3.9) shows immediately that the [a] matrix becomes 
diagonal. Likewise then Eq. (3. 12 ) becomes an explicit finite difference 
operator; The GIM formulation thus allows either explicit or implicit oper- 
ators witli tile finite difference coefficients containing all of the geometric 
transformations. 


The GIM/STAR code contains only the explicit finite-difference analogs. 
A fully implicit scheme for three-dimensional problems requires the 
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simultaneous solution of a large system of algebraic equations. These 
operations are not as readily vectorizable on a machine such as STAR as 
are the fully explicit schemes. Reference 1 contains examples of ♦‘he types 
of explicit schemes which the GIM code can generate. The interpolant 
functions, I, described in this section are general in the formulation. The 
current GIM/STAR code contains only the trilinear inter polants with their 
corresponding orthogonal weight functions. Figure 3-5 shows these inter- 
polant and weight functions that are presently coded in the STAR program. 

The are constants, read into the program, which determine the specific 
type of finite difference algorithm (see Section 4). 

3.2 CODE LAYOUT ON STAR 

The GIM code is split into three separate modules for execution on 
the STAR- 100 system. These are the geometry module, GEOM, the nu- 
merical integration module, IN TEG and the graphics module, GIMPLT. 

Figure 3-6 is a block diagram of the modular structure. All program files 
reside on the STAR Access Station permanent file system in UPDATE format. 

The GEOM and INTEG modules have variable dimensions in the source 
code which are replaced by the actual dimensions by preprocessing the COMPILE 
file before compilation. This is done on the "Z” side by a "dynamic dimension" 
routine. The names of these routines are DYNMAT and DYNDIM for the GEOM 
and INTEG modules, respectively. The input data required for these programs 
is described in the sections on the GEOM module and INTEG module. DYNMAT 
and DYNDIM read the compile file on unit 8 and write the "processed" compile 
file on unit 3. Therefore the UPDATE command must have the parameter 
C = TAPES specified. TAPE 3 is then passed over to STAR via the TOSTAR 
command for compilation. 

The GEOM and DMTEG modules are compiled and executed on STAR, 
hence the "processed" compile files for these modules must be passed over 
to STAR via the TOSTAR command. The GIMPLT module is executed on 
the "Z" machine itself so that no TOSTAR command is needed. 
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Fig. 3-5 - General Interpolants and Weight Functions Used in the GIM/ 
STAR Code 







Fig. 3-6 - GIM Code Modular Structure on STAR 
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3.3 EXECUTING GIM ON STAR 

Executing GIM on the STAR computer system, involves the preparation 
of several items other than the input data. These include: (1) the dynamic 
dimension routine data; (2) preparation of the LOAD card; (3) use of CREATE 
cards for data files; and (4) the remainder of the runstream which is generally 
standard. The description of the dynamic dimension data is given in the 
sections describing the individual modules. Example runstreams are pre> 
sented and discussed in Section 3.4 which provides a guide in runstream 
preparation for GIM runs. This section deals in some detail with items (2) 
and (3) above. The LOAD card involves the specification of large or small 
pages, page boundary allocation for common blocks and the size allocation 
for the controllee file. These are discussed in detail along with a descrip- 
tion of the data files and how to calculate their size requirements. This will 
determine whether a CREATE card is needed. 

File space is created automatically for small files. Larger files must 
be created by a CREATE card. The current default file allocation at NASA- 
Langley is 128 small pages (512 words each). The file allocation will be 
extended up to three times, adding 85 small pages each time, if the original 
allocation is exceeded. Files which require more than the maximum of 384 
small pages (196608 words) must be created with a CREATE command. (See 
STAR Programming Manual and example runstreams.) 

3.3.1 File Descriptions 

File l6: Work file used in GEOM module only. Contains geometry 
information blocked with 500 nodes (NN) or less, per block. 

Size = 10 * NN 1 words minimum 

If less than 500 nodes per block are used, the total file size will be 
larger since records of length 5001 are used even if they are not full. The 
number of nodes per block is determined automatically by the GEOM module 
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based on the zone/region configuration and Is not controlled directly by the 
user. 

File 17; This file contains the nodal analog coefficients. It is gen- 
erated in the GEOM module and used by the INTEG module. It consists of 
a single unformatted binary record of length (NWDS) given below. 

NWDS = 16*NN + 18+NSP + 6 (2D) 

= 48*NN + 18*NSP + 6 (3D) 

where 

NN is the total number of nodes, and 

NSP is the number of special boundary terms (see GEOM 
and INTEG input guides). 

The record length is adjusted upward to be a multiple of large pages (65536 
words). File 17 is mapped onto the nodal ai^^log common block /Q3MAP/ 
in INTEG via the Q30PNMAP call. (See STAR Programming Manual, Section 
5.2.2.) 

File 20; This file is created in the GEOM module and contains all 
geometry data. It is read by both the INTEG and GIMPLT modules. It is 
a formatted file written using a 6E22.14 format. The number of words (NWDS) 
contained in the file is given by the formula below. 

NWDS = 6*NN (2D) 

= 9*NN (3D) 

Since the file is formatted with 22-character word lengths, the file size is 
increased by a factor of approximately 2.2. 

File 22; This file created by INTfcG contains the flowfield solution 
data.^ It is used as a restart file by INTEG and as a flowfield data file by 
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GIMPLT for plotting velocity vectors and contours. It consists of one or 
more records of length NWDS given below. 

NWDS = 5*NN (2D) 

= 6*NN (3D) 

For a two -gas case, additional NN words are contained in each record. The 
number of records is determined by the iteration increment at which the file 
is written. This file is also formatted (6F22.14), so its length is increased 
by a factor of approximately 2.2. 

3.3.2 Controllee File Sizes 

An adequate length for the controllee file must be specified on the LOAD 
card. This requires that some estimate of the size of the controllee file must 
be made prior to execution. Approximate formulas are given for estimating 
the sizes for GEOM and INTEG. 


GEOM 


Number of 
large pages 


125 * NN j. 1 
65536 

307 * NN X “4 
65536 


(2D) 

(3D) 


The value on the LOAD card is specified in small pages (512 words). The 
GRLPALL parameter is always used in the LOAD card for the GEOM module, 
however, indicating the use of large pages (65536 words, or 128 small pages). 
Therefore, the number of large pages required should be calculated based on 
the number of nodes NN and then multiplied by 128 to obtain the value to use 
on the LOAD card. (See example runstreams in Section 3 4). 


INTEG 


Most of the space required by the INTEG controllee file is taken up 
by seven common blocks. The size, L, of these common blocks is given by 
the set of formulas on the following page. 


3-19 





IIIMIIIB 


Where 


/PRIM/ 

/EBUF/ 

/UBUF/ 

/XBUF/ 

/TAU/ 

/AXSYM/ 

/VTEMP/ 


L. = 6 MN + 
L s: MEQ (2 
L. = MEQ (2 
L = 6 MN + 
L = 7 MN + 
L = MNA (4 
L. = 4 MN + 


MNZ -I- MNS 
MN + MNZ) 
MN + MNZ) 
3 MNZ 
5 MNZ 
+ MEQ) 

2 MNS. 


MN 

= 

total 

number of nodes 

MNS 


MN 

for two gas case 


= 

1 

for 

single gas case 

MNZ 

= 

MN 

if 

3D 


s 

1 

if 

2D 

MNA 

= 

MN 

if 

axisymmetric 



1 

if 

not axisymmetric 

MEQ 


4 

2D 

-single gas 


- 

5 

2D 

-two gas or 3D -single gas 



6 

3D 

-two gas 


These common blocks should be assigned to large pages via the GRLiP 
parameter In the LOAD card. Unless the total space for these common 
blocks occupies less than two or three pages, these blocks should be 
assigned to page boundaries as in Runstream Example 3. For a small 
problem they may be grouped together as in Runstream Example 2. 

The remaining large common block in INTEG, /Q3MAP/, is assigned 
to a large page boundary via the GROL parameter in the LOAD card which 
causes no space to be allocated in the cpntrollee file. File 17, which con- 
tains the analog coefficients, is mapped onto this common block via a call 
to Q30PNMAP in INTEG. 
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The space needed for the INTEG controllee file is then the total number 
of large pages required for the set of common blocks above plus one additional 
large page for the remainder of the program. The space parameter is given 
in small pages, therefore, the number of large pages should be multiplied by 
128. If in doubt, add an additional large page. Any unused space allocated 
for the controllee file will be released. 

3.4 EXAMPLE RUNSTREAMS 

3.4.1 Example Runstream 1 (Fig. 3-7) 

This runstream illustrates execution of the GEOM module by itself for 
a large three-dimensional problem (3D shear layer — 7904 nodes). The output 
files, FILE17 and FILE20, are saved back on the "Z" side via the TOAS com- 
mand. Note that FILE17 is a binary file (FILE17=BI) and F1LE20 is a for- 
matted file. Any updates to the GEOM module would follow the *ID GEOMODS 
card. 

Calculation of File Sizes; 

File 16 

L = (10) (7904) + 1 

= 79041 words minimum 
=: 155 small pages minimum 

use 200 in CREATE card 

File 17 

L = (48) (7904) + (18) (7904) + 6 
= 521665 words 
= 1019 small pages 
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n CM60000 * T I OO • 

UStK,6B0077C. 

CHAKGE* 1o1S^37«LWC# 

GET ( ULDPL-GEUiViAT ) 

GET( OYNMaTsDYNMAT ) 

UPDATE< F,C = TAPc:8 ) 

OYNMAT* 

TOSTAR( Input *TAPE3) 

7/8/9 

*IU GEOMODS 
7/8/9 
7904 3 

7/8/9 

stoke 680077 400SDS GLJiviDECK B 
STKGEOM*t200* 

FOKTKAN( I =TAPE3«hJ = Gt0MB*O = L.B) 

CREATE (File 16^ 200 * T = P ) 

CREATE (File 17* 10 1p»T = P) 

CREATE (F 1 LE20006, T=P ) 

load ( GEO, viB»CN = Gf OmGO » 5 1 20 « GRLPALL= ) 
GEOl'-i GO- 
TO AS ( 2=680077C* F IlE17 = BI -FTlEPO) 

**♦ GEOM DATA *** 

• 6/7/S/9/ 


Fig. 3-7 


GEOM Module Execution for a Large Three-Dimensional 
Problem 


3-22 





File 20 


L 


(9) (7904) (2.2) 
156500 words 
306 small pages 


Controllee File Size Calculation; 

L = (307) (7904) + 3 large pages 
= 2426528 words + 3 large pages 
- 40 large pages 
= 5120 small pages 

Note that the CRLPALL parameter is used in the LOAD card. 

3.4.2 Example Runstream 2 (Fig. 3-8) 

This runstream illustrates the execution of GEOM and INTEG together 

for a small problem (340 node 3D). No files are saved on this run. All file 

sizes are small, thus no CREATE cards are needed. 

Calculation of Controllee File Sizes; 

GEOM 

L = (307) (340) + 3 large pages 
= 104380 words + 3 large pages 

= 5 large pages 
= 640 small pages 

INTEG 

This problem is three-dimensional, single gas. Hence, 


U- 
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GEUI NT*CM6000Ci * T 1 OO • 

UStk*680o77C« 

CHARGE « io 1857,i_WC. 

GET( OLDPl = GEO"/iAT ) 

GET ( DYNMaT = DYNiv]AT ) 

UPDATE(F,C=TAP£8 ) 

DYNMAT. 

COPYCF { TAPE3 ♦ 'jEOMC ) 

REWIND<GeOMC)‘ 

RETURNt OlDPL ) 

RETU«N( TAPE3) 

RETURN! TaPEB ) 

GET! CLDPl= INTEG) 

GET ! DYND I M = DYND 1 ) 

UPDATE ! F , C=TAPE8 ) 

OYND IM* 

COPYCF! TAPE3* INTEGX > 

REWIND! INTEGX) 

RETURN! OlDPL ) 

RETURN! TAPE3 ) 

RETURN! TAPEB) 

TOSTAR! InPUT»GEOMC» iNTiGX) 

7/8/9 

*ID GEOMODS 
7/6/9 

340 3 

7/8/9 

*ID INTMODS 
7/8/9 

340 3 0 340 

7/6/9 

STORE 680077 400SDS TtSTUECK 8 
STRSIDEtJlOO. 

FOKTRAN! I =GEUiviC ♦ o ^ G c. O H B • U = LB ) 

LOAD! GlOM8*CN=GEOmGU*640»GRLPALL= ) 
GEOHGO* 

FORTRAN! I s 1 NTEGX » B= I NT EGb • 0=LB ) 
l_OAD! INTeGB»CN= I NtEGO *256 

♦ GRLP = *PRlM**EbUF ,*UB‘JF *^^XbUF**TAUi*VTEFiP 
*GR0L=*U3MAP) 

INTEGO. 

7/8/9 

*** GEOM data *** 

7/8/9 


*** INTEg DATA *** 


6/7/S/9 


Fig, 3-8 - GEOM and INTEG Module b for a Small Three-Dimensional Problem 
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MN = 340 
MNZ = 340 
MNS = 1 
MNA = 1 
MEQ = 5 


The common block sizes are: 


/PRIM/ 

/EBUF/ 

/UBUF/ 

/XBUF/ 

/TAU/ 

/AXSYM/ 

/VTEMP/ 


L = (6) (340) + 340 +1 
= 2381 words 
L = (5) [(2) (340) + 340] 
= 5100 words 
L = (5) [(2) (340) + 340] 
= 5100 words 
L = (6) (340) + (3) (340) 
= 3060 words 
L = (7) (340) + (5) (340) 
= 4080 words 
L = (1) (4 + 5) 

= 9 words 

L = (4) (340) + (2) (1) 

= 1362 words 


The total for common is 21092 words which is less than one large page. 
Therefore, all six common blocks are grouped on one large page. (COMMON/ 
AXSYM/ is not used.) The controllee file size requirement is two large pages, 
hence 256 small pages. 

3.4.3 Example Runstream 3 (Fig. 3-9) 

This runstream illustrates execution of INTEG for a large problem 
(3D shear layer, 7904 nodes) utilizing files created and saved in a previous 
GEOM execution. Note that FILE 17 and FILE 20 are attached as F1LE17A 
and FILE20A. When files are saved on the access station ("Z" side) a letter 
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lNTc.Gi'''*Cf'^600CJ»T^"'0* 

JSl * 6bCr7"7C • 

I 1 057 * l_r<C • 

GET ( CUDPL- irvTEG) 

GET ( u YND I iVi = DYND 1 M ) 

UPDATE (F,C=TAP£0) 

DYND I 

ATTACH(FILE17 = F1i_e17A) 

ATTACH(F1LE20=FILE20A) 

TOgTaP( InPUT%TAPE 3*F I Lt 17=B I //U * F I Lt20 ) 

7/0/9/ 

-5<-ID MCO^ 

7/8/9 

STUPE 680077 ^OOSOS INfEGDCK B 
::jTk’lNTE»TL50« 

FCPTPAN( I -TAPE3»8 = IiMTtGB*0=LB) 

LOAU( INTEGb*CN= 1 NtEGO ♦ 1536 

^ GklP = *PR I ivi ,Gi->.L.P= l£b JH »GKLP = *UBUF * GPi-P = *XbUF *GPLP = *TAU« GRl_P = *VTEMP 
» GPUL =*U3MAP ) 

I NTtiGU. 

TOA0 ( Z=6S0077C*F IlE22 ) 

7/0/9 

*** DATA *** 


6/7/S/9 




is added to make the file name unique if a file by that name already exists. 
Note also that FILiE22, the restart file, is being saved via the TOAS com- 
mand. Following are the common block length calculations and controllee 
file size calculation. 


/PRIM/ 

L 

/EBUF/ 

L 

/UBUF/ 

L 

/XBUF/ 

L 

/TAU/ 

L 

/VTEMP/ 

L 


(6) (7904) + (7904) + 1 
55329 words ^0.85 large pages 

(5) [(2) (7094) + 7904)] 

118560 words ^1.81 large pages 
same as /EBUF/ 

118560 words J^l. 81 large pages 

(6) (7904) + 3 (7904) 

71136 words ~1.09 large pages 

(7) (7904) + 5 (7904) 

94848 words ^^1.45 large pages 
(4) (7904) + 2 

31618 words %0.48 large pages 


Since the total requirement for the common storage is approximately 7.5 
large pages, these common blocks are assigned to large page boundaries 
via the GRLP parameters. Since die common blocks are assigned to page 
boundaries, the total space required is 10 large pages for the common plus 
1 for the remainder of die program, or 11 large pages. This converts to 
1408 small pages. On a problem this size, allow 1 extra large page in the 
controllee file to be sure it is large enough. Hence, use 1536 as the length 
in the LOAD card. 

3.4.4 Example Runstream 4 (Fig. 3-10) 

This runstream illustrates the execution of the GIMPLT module utilizing 
files created and saved on a previous execution. Note that the GIMPLT module 
resides entirely and is executed on the STAR Access Station. Nd TOSTAR 
command is needed. The TAPE22 attach statement is not needed if grid 
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tiLDG12A7A Jl_ HUNT 


PLOTXX«CM 1 1 1000* Tl 00 • 
USEk,660077C. 

CHAkGE. 1018S7.LWC. 

GET ( OLDPl=PLOT I ) 

UPDATEIF) 

KTN( I =C0MP lLt£tL = 0 ) 

ATTACH( TAPE20=F1LE20 ) 

ATTACH( TAPE22=FILe22 ) 
ATTACH(LRCG0SF/UN=L1BRAWY ) 

«FL (111 000 ) 

LDScT (LlB^LkCGOSF , P K £ St T = I NUEF ) 
LGO. 

ATTACH( PlOT/UN^L IBRARY ) 

PLOT. VAkr AN 
7/S/9 

*ID KORChG 
*-D ^ORE.i 

common A(210AO) 

*D K0RE*2 

KMaX=2 1 040 

7/8/9 

*** DATA *** 

6/7/S/9 


Fig. 3-10 - GIMPLT Example Run Stream 
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plots only are being generated. This execution illustrates the use of the 
Varian post processor. Other post processors are available and the user 
should refer to the Graphic Output System User^s Guide for details. 

The example illustrated is for a model with 2000 nodes. The calcu- 
lations necessary for the field length (memory requirements) and the com- 
mon block dimension are given below. Refer to Section 6.5 for discussion 
of these calculations. 


NX 


2000 nodes 

KMAXjq 

= 

5040jq + 8*NX 


= 

5040jq + (8) (2000) 



21040^0 

KFLio 

= 

KMAX^q + 16000^q 


= 

21040 Q + 16000 Q 


= 

37040,„ 

Choose CM 


lllOOOg 


= 

37376jq > 37040^q . 


The A array dimension in blank common, KMAX, is calculated to be 
21040. The runstream illustrates the UPDATE cards to set this dimension. 
The' corresponding field length requirement is lllOOOg. This must be set 
in two places: the job card CM parameter, and in the RFL card immediately 
preceding the LDSET card. 
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4. GEOMETRY MODULE GU IDE (GEOM) 

4. 1 USE OF THE MODULE 

Module 1 of the GIM/STAR code calculates the geometry, grid and 
nodal analogs for the problem. This module is always executed first and 
produces two output files of data for use in Modules 2 and 3. The GIM code 
GEOM module as implemented on the STAR system at NASA -Langley has a 
'•dynamic dimension" capability. A preprocessor, called DYNMAT, must be 
executed with GEOM to provide dimensions to the problem. Figure 4-1 
illustrates the logic sequence for using Module 1. 

The output files from GEOM are the following; 

• File 20 — Contains the grid point description 

• File 17 — Contains the nodal analog matrices 

If GEOM is executed in a runs tr earn, in a stand alone mode, these files must 
be saved for later input to the INTEG and/or GIMPLT modules (see Section 3). 

The execution time for the GEOM module is very strongly dependent 
upon the dimensionality of the problem and the number of nodes used. The 
following information is provided from execution of the two example problems. 


Configuration 

Nodes 

User 
CP Sec. 

CRU 

2-D Scramjet 

940 

27 

22 

Simulation 




3-D Scramjet 

7904 

1009 

449 

Simulation 
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GEOM Module 
















Thus for small, two-dimensional problems, the execution time for GEOM is 
negligible comt>ared to INTEG. However, for larger, three-dimensional 
cases, the GEOM run time can be comparable to the INTEG run times. 

Further experience with the STAR version of GIM should provide more 
information on run times. 

The following subsections describe the input data and output listings 
for this module. Section 4.2 describes the input for the DYNMAT deck and 
Section 4.3 gives a "mini-manual" of the input format for GEOM. A de- 
tailed discussion of the input variables is provided in Section 4.4. The 
output listings are described in Section 4.5. 

4.2 DYNMAT DECK 

The GEOM module of the GIM/STAR code has a "dynamic -dimensioning" 
feature. The larger dimensional variables and arrays have parameters in 
place of actual integers. The DYNMAT program uses the UPDATE feature 
of the CDC software to replace these parameters with real dimensions so that 
GEOM can be compiled. These dimensions are under user control via input 
to the DYNMAT program. Execution of this "pre-processor" routine must 
be made before executing GEOM as explained in Section 3.4. This section 
describes the input data for the DYNMAT program. 

The required data are input on one card with the following format: 

FORMAT (215) . 


NX 

The total number of nodes in the problem. All flow variables are 
dimensioned by NX. 
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IDIM 


The dimensionality of the problem. 

IDIM = 2 for two-dimensional geometries 

= 3 for three-dimensional geometries. 

This flag controls the number of coefficient matrices which are dimensioned. 

The remaining dimensions that are required are calculated internal 
to DYNMAT from these two input parameters. 

4.3 INPUT CARD SUMMARY 

This subsection presents a summary of the input cards and formats 
for the GEOM module. A description of each input parameter and its options 
aredetailed in the following subsection (4.4). After a user becomes familiar 
with the GEOM inputs, this summarized input gqide can be used to quickly 
identify each card and its contents. 

Three basic formats are used to input the data to GEOM: 

ALPHANUMERIC A 6 

INTEGER 15 

DECIMAL E10.4 

Integers are dius right justified in five column increments. Decimal, or 
floating point data occupy ten columns each with, preferably, a decimal point 
punched on the card (see examples. Section 8). 

Card Type Parameter List/Format 

1 HEADER(I), 1= 1,72 

(12A6) 

2 NZONES, IDIM, ISTEP. IMATRX, IMATE 

(515) 
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Card Type 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 


Parameter LLst/Format 

IWRITE, LWRITE, NWRITE 
(315) 

ALPHA(I). 1= 1.4 (or 8) 

(8E10.4) 

NSECTS 

(15) 

MAPE(I), 1= 1, 12 
(1215) 

MAPS(I), 1= 1, 6 
(515) 

(IBWL(I), 1= 1,6), ITRAIN 
(715) 

(NNOD(I), 1= 1, 3), (ISTRCH(I), 1=1, 3) 

(615) 

DIVPI(I), 1= 1, 3 
(3E10.4) 

[AETA(J,I), 1= 1. NNOD(J)] , J = 1, IDIM 
(8E10.4) 

[(AC(I,K, J), 1= 1,8), J= 1,4 or 12], K= 1,5 
(8E10.4) 

[AS(I, J), 1= 1,8] , J = 1,6 
(8E10.4) 

(PT(I, J), 1= 1, 5), J= 1, 4 or 12 
(8E10.4) 

[(PMAX(I,K, J), 1= 1, 5), ETAMAX(K, J), K= 1,4] , 

J = 1, 4 or 12 

(6E10.4) 


Notes; Card types 14, 15 are repeated for each point and edge. 

Card types 6 through 15 are repeated for each of NSECT sections 
in a zone 

Card types 12, 13, 14, 15 are not re -input for points, edges and 
surfaces in common with the preceding section. 

Card types 5 through 15 are repeated for each zone of NZONES 
to be processed. 
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Parameter List/Format 

NDX, NDY, NDZ, ISNOPT 
(415) 

N1.IC,NT 
(315) 

4.4 DESCRIPTION OF INPUT DATA 

This subsection presents a description of the input parameters listed 
in Section 4.3. Each parameter is identified as to its usage in the~C3EOM 
module with options of each shown. Reference to the figures and tables 
must be made to explain some of the input parameters and order of input. 

All of the card types are not necessarily input for a specific case. Much 
of the input depends on the dimensionality of the problem (two- or three- 
dimensional). Also, certain of the input flags on early cards dictate which 
of the latter ones are read. Section 8 gives examples of the GEOM input, 
which can be used as a guide for setting up a new problem. 

In the subsequent description of the input parameters, a nomenclature 
is used for describing the geometry that is key to understanding the sequence 
of events. The key words are 

ZONE, SECTION and SEGMENT. 

The full flow domain may be input as one or more zones. A "zone" is herein 
defined as a first level partitioning of a full domain for the purpose of analyt- 
ically describing a shape. A full domain is thus partitioned into "sub-domains" 
called zones, such that each zone can be readily defined by suitable analytical 
functions. 


Card Type 
16 

17 


A "section" is defined herein as a second level partitioning of a full 
domain or equivalently, a first level partitioning of a zone. This further 
breakdown is done to facilitate easier handling of geometrical changes. 
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A problem, for example, may be best described by one zo 
subdivided into many sections. 

Zones can be set up independently and either left th 
with other zones. Thus, a building block modular concep 
completely general geometries can be treated. 

A "segment" is defined as a partitioning of an edge, 
general, for the purpose of grid distribution. An edge of 
sist of one segment in which case the nodal distribution w 
segmenting an edge into more than one nonuniform piece, 
irregular grid patterns is made easier. 

Three other names of extreme importance in under 
input are: 


POINT, EDGE and SURFACE. 

Figure 4-2 illustrates each of these for a general hexahe 
defined as the x, y, z coordinates of the corner of the genj 
An edge is a line segment connecting two points. A surf; 
shape formed by connecting four points with an analytical 
should be studied before proceeding with the card-by-cai 

CARD TYPE 1 Format (12A6) 

Problem Identification Labels 


HEADER 


Any alphanumeric information that the user wants 
specifics of the problem. Columns 1-72 are read and pi 

CARD TYPE 2 Format (515) 

Problem Option Control Flags 
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Surf 4 



PT 


6 


Note: , Surf 1 — Back Wide Face 
Surf 2 — Bottom 
Surf 3 — Front Wide Face 
Surf 4 — Top 

Surf 5 — Back Narrow Face 
Surf 6 — Front Narrow Face 


(^1* 


0) 

iriy 

0, q 

3^ 



1) 


l,q 

3^ 

(0,r] 

2’^ 

3^ 

(1. 11 


3^ 


Fig. 4-2 - General Hexahedral Showing Point, Edge and Surface 
Numbering Sequence 
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NZONES 

Number of zones into which the full domain geometry is broken. 
Regular regions can usually be input as one zone, but highly irregular shapes 
require NZONES > 1. (No limit.) (See Section 8 for examples of zones, 
sections and segments.) 

IDIM 

Flag indicating the dimensionality of the problem being run. 

IDIM = 2 Two-Dimensional 

= 3 Three-Dimensional 

The INTEG module has an option for IDIM= 1 which is axisymmetric flow. 
The GEOM module does not need to know if the flow is to be planar or 
axisymmetric. Only 2, 3 are used for IDIM in the GEOM module. 

ISTEP 

A flag indicating whether a one-step or two-step time integration 
scheme is to be used. 


ISTEP = 1 One Step 
= 2 Two Step 

The GEOM module will output either one or two sets of finite difference 
matrices depending on ISTEP. 

IMATRX 

A flag which allows only a geometry and grid to be generated or both 
a grid and the finite difference matrices. 


L. 


4-9 



IMATRX = 0 Provide both a grid network 
and matrices 

= 1 Only compute grid points 
and nodal connections; do 
not assemble matrices 

Common use of the GEOM module is to set IMATRX= 0. For some large 
problems, the user may want to wait until a good grid is obtained before 
assembling the matrices. For current use on STAR, it is suggested that 
the zero value be used. 

IMATE 

GEOM has the option of processing more than one zone in a given run. 
The IMATE parameter allows the user to either mate the zones together to 
form a full domain; or to keep the zones separate — unmated. 

IMATE = 1 Mate zones together 

= 0 Do not mate the zones 

Ordinarily, the user will want IMATE = 1 to form a full region. For some 
complex shapes/grids, certain zones may need to be kept separate. (See 
Section 8 for examples of zones, sections and segments.) 

CARD TYPE 3 Format (315) 

Output Print Options 

IWRITE 


A flag used for debugging printout only. 

IWRITE = 0 No debug print 

= 1 Provides intermediate output 

The user should probably set IWRITE = 0 for all cases unless a problem 
occurs during execution. 
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LWRITE 


An option controlling the printout of each element matrix. 

LWRITE = 0 Do not print element matrices 

= N>0 Print element matrices for 
every element 

The user will probably use LWRITE = 0 most of the time unless a problem 
occurs during execution of the GEOM module. 

NWRITE 

This flag is the more important of the three print options. This con- 
trols the amount of grid point printout that is given for a run. 

NWRITE = 0 Print only those grid points 

which lie on a boundary 
th 

= N>0 Causes every N node point 
to be printed out 

The output format of the GEOM data is described in Section 4.5. For large 
problems with many thousands of nodes, NWRITE should be chosen carefully 
to restrict the amount of paper produced on STAR. A plot of the grid is 
usually more instructive than a massive nodal printout. 

CARD TYPE 4 Format (8E10.4) 

Weight Function Coefficients, 

ALPHA(l) 

These coefficients are used to define the weight functions and hence 
the finite difference method. The reader may want to review the theory in 
Section 3 and Ref. 1 for definition of the ALPHA parameters. 
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Table 4-1 gives values for AJLPHA which produce several useful 
schemes, such as MacCormack, etc. 

For two-dimensional cases, only four values of a are input for each 
step of the method. All values are input on one card in Format (8E10.4). 

For one-step methods, only input four values of ALPHA; for two step methods, 
input eight values of ALPHA. 

For three-dimensional problems, eight values of a are needed for 
each step. One step methods require all eight values on one card. Two 
step methods requires two cards type 4 with eight values each. 

CARD TYPE 5 Format (15) 

Number of Sections 


NSECTS 

Number of sections within a zone for the zone currently being input. 
Recall that a flow region may be broken into zones. Further, each zone 
can be subdivided into a number of sections such that analytic functions 
describe the shape of each section. There is no limit to the value of 
NSECTS within a zone. 

CARD TYPE 6 Format (1215) 

Edge Shape Function Indicators 

MAPE(I) 

These are integer flags that define the edge shape functions for the 
section being read-in. The edges are input in nvunerical order according 
to Fig. 4-2. The user should study this figure before setting up an edge 
and surface input deck. 

The values of MAPE depend on the type of analytic function describing 
the edge. The library of edge shape functions in the GIM/STAR code are 
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Table 4-1 


VALUES OF alpha FOR SEVERAL USEFUL FINITE 
DIFFERENCE SCHEMES (CARD TYPE 4) 


Forward-F orward-(F orward){FFF) 


2 

i-D 

3 

-D 

a 

Value 

a 

Value 

1 

1.0 

1 

1.0 

2 

l.OE-8 

2 

l.OE-8 

3 

l.OE-16 

3 

l.OE-12 

4 

l.OE-8 

4 

l.OE-8 



5 

l.OE-8 



6 

l.OE-12 



7 

1. OE- 16 



8 

l.OE-12 


Backward-Backward- (Backward)(BBB) 


2 

.-D 

3 

-D 


Value 

a 

Value 

1 

l.OE-16 

1 

l.OE-16 

2 

l.OE-8 

2 

l.OE-12 

3 

1.0 

3 

i.OE-8 

4 

l.OE-8 

4 

l.OE-12 



5 

l.OE-12 



6 

l.OE-8 



7 

1.0 



8 



l.OE-8 
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Table 4-1 (Con.clu.ded) 


Backward -F orv/ard-(F orward)(BFF) 


2 

:-D 

3- 

D 

a 

Value 

a 

Value 

1 

l.OE-8 

1 

l.OE-8 

2 

1.0 

2 

1.0 

3 

l.OE-8 

3 

l.OE-8 

4 

I.OE-16 

4 

1.0E-X2 



5 

l.OE-12 



6 

l.OE-8 



7 

l.OE-12 



8 

I.OE-16 


I 

Centered Scheme 2-D 

Forward-Backward 2-D 

a 

Value 

a 

Value 

1 

1.0 

1 

l.OE-8 

2 

-1.0 

2 

I.OE-16 

3 

1.0 

3 

I.OE-8 

4 

-1.0 

4 

1.0 


Note: Standard MacCormack scheme will use FFF 
on step 1 and BBB on step 2 , 
























listed in Table 4-2. Three-dimensional problems require all 12 edges while 
two-dimensional problems require only edges 1, 2, 3 and 4. If any shape 
function other than a linear segment is used, then coefficients must be input 
on card type 12 to define this function. 

An option is provided to allow distribution of grid points in each co- 
ordinate direction. This is done by segmenting each edge. Each of the 
edges may consist of up to five segments with each of these segments having 
its own shape f\inction. The values of MAPE then can consist of up to five 
integers packed into one word MAPE (I). The segment indicators are input 
in chronological order of increasing 77 for each edge with the final packed 
integer being right adjusted. For example, if MAPE(4) = 312, then edge 4 
consists of three segments; the first segment being type 3, the second seg- 
ment type 1 and the third segment type 2 , If a single shape function describes 
an edge, the only one indicator is used, right adjusted. The example prob- 
lems of Section 8 show how this segmenting procedure can be used. 

CARD TYPE 7 Format (6l5) 

Surface Shape Function Indicators 


MAPS (I) 

These flags are input only for three-dimensional problems since two- 
dimensional geometries are defined completely by the edge functions. 

The input sequence for the surfaces are shown in Fig. 4-2 for designa- 
tion 1 to 6 , The surfaces must be input in this order. The values for MAPS 
are given in Table 4-3. There are three basic types; flat plate, cylinder and 
edge of revolution. The latter case requires input of coefficients on card 
type 13 to define the revolution axis, etc. The user should refer to Fig. 4-2, 
Table 4-3 and the examples given in Section 8 . 
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Table 4-2 


EDGE SHAPE FUNCTION INDICATORS AND 
COEFFICIENT DEFINITIONS 


Type 

Map 

(Card 6 ) 

Coefficients 
(Card 12) 

Linear 

1 

None. Omit Card 12 

Circular Arc 

2 

ACj^ — AC^ are (x,y,z) coordinates 
of center of arc. AC^ is the ex- 
pansion angle from nozzle, if any. 

Conics (Parabola, 
Ellipse, Hyper- 
bola) 

3 

AC, — AC^ are (x, y,z) coordinates 
of focal point. AC 4 — AC^j are 
(x, y, z) components of unit vector 
along axis of conic. 

Helical Arc 

4 

ACi— AC^ are coordinates of 
center of arc. AC^— AC^ are 
components of unit vector along 
axis. 

T rignometric 
F unction of X 

e 

5 

ACj^ — AC^ arc coordinates of 
reference point corresponding 
to Xg= 0 . AC 4 — AC^ are used in 
equation 

y = Yq + AC^ sin (AC^X^ + AC.^) 
+ ACg cos (AC^X^+ACy) 

T rignometric 
Function of 9 

6 

AC^— AC^ are coordinates of 
center of sweep angle 6 . AC^-*- 

ACg are used in equation 

y = AC 4 + ACg sin (AC.^ 0+ACg) 

+ ACg cos (AC.^ 0 + ACg) 

Special Functions 

7 

The user can modify the GEOM 
code to include any special function. 
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Table 4-3 


SURFACE SHAPE FUNCTION INDICATORS 
AND COEFFICIENT DEFINITIONS 


Type 

Map 
(Card 7) 

Coefficients 
(Card 13) 

Flat Plate 

1 

None. Omit Card 13 for this 
surface. 

Cylindrical 

Z 

None. Omit Card 13. 

Special 

Application 

3 

The user can supply special 
function by modifying the 
code for a specific application. 

Edge of 
Revolution 

4 

Surface formed by revolving edge 
about an axis. ASj^—AS^ is origin 
of surface coordinate system. 

AS^ — AS^ are components of unit 
vector along axis of revolution. 
AS.^ indicates which direction 

the edge will revolve in. 
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CARD TYPE 8 


. Format (715) 


Boundary Indicators and Node Numbering 
Sequence Specification 


IBWL(I) 

These integer indicators are used to tell the INTEG module which 
boundary conditions to apply to an edge/surface. Table 4-4 lists the allow- 
able boundary types. 

IBWL(I) is the flag for edge I or surface I. If IDIM= 2, then the IBWL 
flags are for edges and if ID1M= 3 then IBWL is for surfaces. For two- 
dimensional problems, there will be four edges to define, 1=1,4. For three- 
dimensional problems, there are six surfaces to define, 1= 1,6. 

Each value of IBWL(I) contains a pac ked integer which contains two 
flags. The right-most integer is the boundary indicator listed in Table 4-4. 
The single digit integer just to the left is a flag which designates the zone 
number to which the edge is to be mated. For example, 

IBWL(3) = 14 

says that edge 3 is a boundary type 4 (wall tangency) and that edge 3 is to be 
mated with zone 1 of a previous input sequence. 

The boundary flags are input in numerical order according to the edge/ 
surface description in Fig. 4-2. Columns 4-5 thus describe edge 1, columns 
9 and 10 for edge 2, etc. 

ITRAIN 


This flag is always input in column 35 of card type 8. The value of 
ITRAIN specify the node numbering sequence in the coordinate directions. 
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Table 4-4 


BOUNDARY NODE DESIGNATORS AND THEIR 
MEANING FOR USE IN THE INTEG MODULE 


IBWL 

Boundary Type 

Remarks 

0 

Constant Node 

All flow variables are held fixed 
at the input value. 

1 

Axis Node 

Designates the axis points for 
axisymmetric flow. 

2 

No-Slip/Stag- 
nation Node 

The density times velocity values 
are held fixed, (pu,pv,pw) but 
pressure, energy and density 
equations are integrated. 

3 

Corner Node 
(3-D) 

For internal 3-D corners, the 
flow is forced to go down the 
corner and hence tangent to 
both planes. 

4 

Free Slip/ 
Tangency 

Velocity vectors are forced to 
be tangent to a wall (or = 0). 

5-7 

Not Currently 
U sed 

The user may modify INTEG 
for special cases. 

8 

One-Sided 

Differences 

A boundary is computed from 
only information inside the 
region. Good for supersonic 
downstream condition. 

9 

Interior Node 

No constraints. All flow 
variables computed from 
the equations. 
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The table below shows the values of ITRAIN for numbering along each of 
three coordinate axes for two- and three-dimensional problems. 


IDIM 

ITRAIN 

Numbering Sequence 

2 

1 

r?2»ni 

2 

2 

'>r”2 

3 

1 


3 

2 


3 

3 



An example of the use of this table is: IDIM = 2, ITRAIN = 1; all nodes in 
the r)^ direction are numbered first at the initial line. Then the nodes 
are numbered along 3it the next station etc. 

CARD TYPE 9 Format (6l5) 

Node Point Distribution Indicators 


NNOD(I) 

For 1= 1, 2, 3, NNOD(I) is the number of nodes in the rjj direction for 
the current section being input. The limit is 100 nodes in each coordinate 


direction or 10 nodes total. 
ISTRCH(I) 

An option is provided to allow 
of the coordinate directions. 

ISTRCH = 0 

= -1 

= +1 

= NNOD 


unequal distribution of nodes in each 

uniform spacing 

reduce spacing in rjj direction 

increase spacing in rjj direction 

user input location of nodes 
on card type 1 1 . 


4-20 



r 


For uniform spacing. Cards 10 and 11 are not needed and are not input. 
If ISTRCH = +1» Card 10 is needed but not Card 11. If ISTRCH = NNOD, then 
Card 10 is not needed, but Card 11 is required. 

CARD TYPE 10 Format (3E10.4) 

Stretching Function Values 


DIVPl(I) 

The values for each coordinate direction which are used to determine 
the degree of stretching by the TANGENT transformation given in Fig. 4-3. 
Typical values of 2.0 and 3.0 are shown as to the relative uniformity of a 
grid. Large values of DIVPI produce essentially uniform spacing while a 
small value gives drastic stretching of the grid points. This card is input 
if and only if ISTRCH = + 1 on card type 9. There are either two or three 
values depending on IDIM= 2 or 3. 

CARD TYPE 11 Format (8E10.4) 

Optional Grid Point Input 


AETA(I) 

The actual rj coordinates of the grid for the section being 

constructed. These are input if and only if any of the ISTRCH (I) were read 
in as NNOD(I). For example, the r] ^ coordinate line could be equally spaced, 
the coordinates calculated with a stretching transformation and the rj^ 
coordinate line could be read into AETA(I), 1= It NNOD(3). 

The AETA values are read in the order rj j then then rj^ as required 
in Format 8E10.4 using as many cards as required to get them all read. The 
AETA values are normalized to range from 0.0 to ,1.0 and hence give a rela- 
tive spacing of nodes. The AETA values are not the actual x, y, z coordinates 
For three-dimensional cases, there are up to three sets of AETA for each 
coordinate line. For two-dimensional cases there are up to two sets of 
AETA cards. Remember to input card type 11 if any of the ISTRCH(I) are 
set to NNOD(I) on card type 9. 
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N/NNOD(I) 


Note: Stretching functions are tangent functions of the form: 


• for reducing nodal spacing 




1.0 


tan 

7T NNOD(I) - N V 

DIVPl(I) ' NNOD(I) ' 

tan 

7T .NNOD(I) - 1.0. 

DIVPI(I) ^ NNOD(I) ’ 


N=l,2,...,NNOD(I) 


• for increasing nodal spacing 


tan 

IT 

( )1 

;Divpi(i) 

' NNOD(I) 

tan 

7T 

,NNOD(I) - 1.0 . 

DIVPI(I) 

^ NNOD(I) ^ 


Fig. 4-3 - Stretching Function Illustrated 
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CARD TYPE 12 


Format (8E10.4) 


Coefficients for Edye Shape Functions 

AC 

The ACs are coefficients which describe the shape functions for each 
edge of the section being constructed. If an edge is segmented, the coef- 
ficients for each segment are input on separate cards in the same order as 
the indicators on card type 6. Card type 12 is input if and only if MAPE > 1 
on card type 6. MAPE = 1 defines the edge as linear and no coefficients are 
required. 

Table 4-2 defines the AC values to be input for each type of edge. The 
examples of Section 8 should clarify the AC input sequence. 

CARD TYPE 13 Format (8E10.4) 

Coefficients for Surface Shape Functions 
AS 


The AS values are the coefficients which define the shape functions for 
each surface of the section being constructed. This card may be input if and 
only if 1DIM= 3 since surfaces are not present in two-dimensional geometries. 
Table 4-3 gives a description of the AS coefficients for each type of surface. 

Note that AS coefficients are required only if MAPS > 2 in Table 4-2. 

For linear or cylindrical surfaces^ no coefficients are necessary to describe 
the surface. Only for edges of revolution (MAPS = 4) or special cases (MAPS = 3) 
is a card type 13 required. 

Each surface which has MAPS > 2 on Card 7 is input on a separate 
card in the same order as the surface indicators on card type 7. 
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CARD TYPE 14 Format (5E10.4) 


CoordLnates of Points 


PT(I« J) , J = 1, Number of Points 
1= 1. 5 

The PT values are the coordinate points which define each corner of 
a general hexahedral. Figure 4-2 shows this configuration with the points 
numbered from 1 to 8. 


PT(1, J) 
PT(2, J) 
PT(3, J) 
PT(4, J) 

PT(5, J) 


the X coordinate of point J 

the y coordinate of point J 

the z coordinate of point J 

flow angle in the x-y plane 
at point J 

flow angle in the x-z plane 
at point J 


For two-dimensional problems there will be four cards of type 14 to 
be input. For three-dimensional problems there are eight cards of type 14 
to be input. 


Important Note; All cards type 14 are not input consecutively. They 
are grouped with cards type 15. See Table 4-5 for 
the exact sequence of types 14 and 15 and the order 
in which they must be input. Section 8 has examples 
of this usage. 

CARO TYPE 15 Format (6E10.4^ 

Segment Extremals for Edges 

PMAX(1,K, J) 


The extremal x coordinate, in the r) ^ 
of edge J. 


direction, for the K 


th 


segment 
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Table 4-5 


INPUT SEQUENCE FOR CARDS TYPE 14 AND 15 
TO DESCRIBE ALL POINTS AND EXTREMALS 
OF SEGMENTS 


Order 

Card Type 

Description 

1 

14 

Point 1 

2 

15 

Extremals for Edge 1 

3 

14 

Point 2 

4 

15 

Extremals for Edge 2 

5 

14 

Point 3 

6 

15 

Extremals for Edge 3 

7 

14 

Point 4 

8 

15 

Extremals for Edge 4 

9 

15 

Extremals for Edge 5 

10 

14 

Point 5 

11 

15 

Extremals for Edge 6 

12 

14 

Point 6 

13 

15 

Extremals for Edge 7 

14 

14 

Point 7 

15 

15 

Extremals for Edge 8 

16 

14 

Point 8 

17 

15 

Extremals for Edge 9 

18 

15 

Extremals for Edge 10 

19 

15 

Extremals for Edge 11 

20 

15 

Extremals for Edge 12 


Note: For two-dimensional geometries, the order 

of input stops after 8 since no other points 
or edges exist. 
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PMAX(2, K, J) 


The extremal y coordinate, in the rj^ direction, for the k“^ segment 
of edge J. 

PMAX(3,K, Jj 

th 

The extremal z coordinate, in the rj^ direction, for the segment 

of edge J. 

PMAX(4,K, J) 

th 

The flow angle in the x-y plane at the extremal point on the seg- 
ment of edge J. 

PMAX(5,K, J) 

th 

The flow angle in the x-z plane at the extremal point on the K seg- 
ment of edge J. 

ETAMAX(K, J) 

th 

The maximum value of the rj coordinate on the K segment of edge J. 
Notes on Card Type 15; 

Each edge may be segmented up to five times. Therefore, cards type 
15 are repeated for each successive segment on edge J. Each segment must 
be input on a separate card type 15. The extremal for the final segment of 
an edge is not to be input since this point is already defined by the PT input. 
The number of cards type 15 for each edge will thus be one less than the 
number of segments on that edge. In particular, if an edge consists of only 
one segment, iw cards of type 15 are input for that edge. 

Table 4-5 should be consulted for the input order of cards type 14 and 
15. For two-dimensional problems, die order of input stops after 8; for 
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three-dimensional problems the order goes up to 20 as shown in the table. 
The points are input on a single card, followed by up to five cards containing 
the extremals. Be sure to note the change of order of the input pairs after 
order 8. Up to order 8, the points are input first — then the extremals. 

After order 8, the extremals are input first — then the points. At order 17> 
the input of card type 14 ceases since all points, have been input. Only ex- 
tremals are input for orders 17-20. 


Notes on Input Sequence of Cards Type 5 Through 15 

• Cards type 6-15 are repeated for each of NSECT sections in a 
zone. 

• Cards type 5-15 are repeated for each zone of NZONES to be 
constructed. 

• Cards type 12-15 are not re-input for points, edges and surfaces 
in common with the preceding section. Since these have already 
been input, the code transfers the entries from one section to 
another. 

• The first section of each new zone however, must be input since 
zones are considered to be independent. Zones can of course be 
mated as previously discussed. 

• When inputting card type 10 for successive sections, only the 
value of DIVPI (ITRAIN) may be changed among the sections. 

The other DIVPI must be the same as for the preceding section. 

• When inputting card type 11 for successive sections, only the 
value of AETA (ITRAIN) may be changed among the sections. 

The other AETA must be the same as for the preceding section. 


CARD TYPE 16 Format (415) 

Coefficient Matrix Parameters 

NDX 

This parameter is the nodal decrement in the x coordinate direction 
for the full flow domain. NDX is used by the STAR/GEOM module to per- 
form the matrix diagonalization. NDX is thus the difference in node numbei 
of adjacent nodes along the x axis. See the two-dimensional Scramjet simu- 
lation case in Section 8 for an example. 
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NDY 


The nodal decrement in the y coordinate direction for the matrix 
diagonalization. 

Nbz 

The nodal decrement in the z coordinate direction for the matrix 
diagonalization. Set NDZ = 0 for two > dimensional problems. 

ISNOPT 


The parameter is used to provide an option on the treatment of off« 
diagonal terms in the matrix handling on STAR. If 

ISNOPT = 0 

then the GEOM module will set the number of special boundary terms equal 
to the number of nodes. In general this will produce a dimension greater 
than is needed. This option should be used only if the INTEG module is to 
be run Immediately following the GEOM module. If 

ISNOPT = 1 

then the GEOM module will calculate the exact number of special off-diagonal 
terms and print out this number. The number can then be input to the INTEG 
module to produce an exact dimension for the special terms. This option 
should be used if the GEOM module is to run in a separate execution stream 
before running the INTEG module. 

The reason for having this option is to reduce the amount of storage 
in the INTEG module. For small two-dimensional problems, it is unimpor- 
tant, but for larger three-dimensional problems the large page faults can be 
reduced by using ISNOPT = 1. 
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CARD TYPE 17 Format (315) 

Nodal Analog Print Control 


N1 


This is the first node of a sequence for which the nodal analog output 
is desired. 


IC 


The nodal analog print increment to use, starting with node Nl« to 
printout the nodal analogs. 


NT 


The total number of nodes, in this sequence, for which analog printout 
is wanted. 


Notes; 


• Any number of cards type 17 can be input. The last card must 
have a - 1 in columns 4 and 5 to tell the code that no more cards 
are coming, 

• The analogs can be printed out in any order. 

• Example of cards type 17: 

15 1 2 

684 2 5 

-1 

This sequence will produce analog printout for nodes 15, l6, 684, 

686 , 688 , 690 , 692 . 

• Any number of nodal analogs may be printed up to the total num- 
ber of nodes in the problem. 
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4.5 OUTPT’'" DESCRIPTION 

The output of the GEOM module consists of two types: 

• Printed grid and nodal analog descriptions 

• Two save files on which the grid and nodal analogs are stored. 

The output save files from the GEOM module are 17 and 20. File 17 con- 
tains the nodal analogs and File 20 contains the grid points and descriptions. 
These files must be saved by user for later (or immediate) input to the 
INTEG (17 and 20) or GIMPLT (20 only) modules. 

The printed output from GEOM consists of two basic types: 

• Geometry/grid description 

• Nodal analog definitions 

Each of these types is now briefly described. The reader should refer to 
the corresponding figures while reading this section. 

The first few pages of output from the GEOM module consists of sum- 
maries of the input data. The control parameters, the zone and section data, 
edge and surface designators are printed (see Figs. 4-4 and 4-5). 

The major portion of the geometry output consists of node point loca- 
tions and geometric data. The pages are headed by the statement 

MESH POINTS FOR RECORD N . 

Figures 4-6 and 4-7 are examples of this output for the two-dimensional 
Sc ramjet problem. As these figures show, the first column is the node 
number followed by the following list: 

X, Y, Z — Cartesian coordinates of the node 

THXY, THXZ — Flow angle descriptors in the x-y and x-z 
planes for the node 
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Fig. 4-4 - GEOM Input Card Printout/Summary 
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Fig. 4-5 - GEOM Input Card Summary 
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Fig. 4-6 - Example Grid Output (Page 1) 
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Fig. 4-7 - Example Grid Output (Page 2) 
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VOL — A volume parameter associated with the node (average 
only and not a true volume). 

NX, NY,NZ — Geometric parameters which are non-zero only 
for boundary nodes. If BC = 4 , these are unit normal vector 
components to the surface. If BC = 3, these are unit tangent 
vectors. 

BC — Boundary condition types (see Table 4-4). 

The user can control the number of nodes printed on the input cards 
previously described. 

Following the grid point printout is a summary of the matrix diagonall- 
zation data and nodal analog print flags. The actual nodal analogs, the finite 
di^erence coefficient matrices, are printed next. Figures 4-8 through 4-11 
are examples from the two-dimensional Scramjet problem. The analog 
printout is headed by statements of the form: 

NODAL ANALOG FOR X DIRECTION, STEP 1. 

The first integer on a line is the node number for which the analog is being 
printed. Following this node number is a list of the connecting nodes and 
the finite difference coefficient. Refer to Fig. 4-8 (line 1 node 129) while 
the following example is given. 

Example: 
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Fig. 4-8 - Matrix Coefficient — Nodal Analog Output (X, Step 1) 
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Analog Output (Y,Step 1) 
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Fig. 4-10 - Matrix Coefficient — Nodal Analog Output (X,Step 2) 
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Reading from Fig, 4-8 on line 1, we see that the x-derivative at node 129 is 
given by; 


^1 
ax li29 


(-10.895807) Ej 29 + (10.085480) 
(9.8651875) Ej^^ + (0.022070995) E^^^ 


This is a type of forward-forward operator in the finite difference sense. 


Nodal analogs for the y and z coordinates are printed next for step 1. 
These are then followed by the analogs for step 2, if a two-step method is 
being used. (See the remaining Figs. 4-9 through 4-11.) This concludes the 
printout of the GEOM module. 
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5. INTEGRATION MODULE GUIDE (INTEG) 


5.1 USE OF THE MODULE 

Module 2 of the GIM/STAR code performs the actual integration of the 
differential equations. The "derivative taker" matrices from module 1 pro- 
vide the type of finite difference scheme such that the integration module can 
be considered a general purpose equation solver. Boundary conditions are 
treated in this module via input flags specified by the user in the GEOM 
module. The output is a tape containing the flow field at user specified 
time increments. This tape can then be printed, or contour maps obtained. 

The GIM numerical analog has the form: 

[■^mnl [®mnl ^n [^mnl [^mn 1 ^n " ® 

where the matrices [A^n] , etc., are obtained from the output tape of module 
1, It is these coefficients which determine the finite difference scheme to be 
used. Schemes which have been tried to date include Galerkin, MacCormack, 
Hopscotch algorithms, centered schemes and Euler methods. 

The calculation can start at time =: 0 with initial conditions as input or 
it can be restarted from the output solution of previous iterations. The 
solution proceeds for a specified number of time integration steps. 

Module 2, INTEG, must be executed after the GEOM, Module 1, has 
been run. Also, the DYNDIM deck must also be run on the 6000 side of the 
STAR system to provide dimension information to INTEG. Figure 5-1 shows 
the logic path of execution for the INTEG module. The requirements are 
Files 17 and 20 from the GEOM Module 1 and the dynamic dimensioning 
infornaation from DYNDIM (discussed in the next subsection). 
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DYNDIM 

6000 



Fig. 5-1 - Logic Sequence for INTEG Module 
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After thLs "preprocessing" is done, the INTEG module is executed on 
STAR and reads data cards in subprogram INPUT. For a restart case. 

File 22 from a previous execution is also input. The code is then executed 
in a loop as shown in Fig. 5-1 until the maximum number of iterations is 
reached. 

This module of the GIM/STAR code usually requires the largest run 
time of the code. The actual CP time is a strong function of the problem 
itself, the grid construction and the options chosen. The following numbers 
are a rough guideline for execution of INTEG on STAR. The numbers are 
the CP seconds required for one iteration of one node point. Total CP is 
obtained by multiplying these estimates by the number of nodes in the prob- 
lem times the number of iterations being run. 


Mode 

Invisc id 

V iscous 

2- D 

3- D 

0.6E-4 

l.OE-4 

l.OE-4 

1.8E-4 


The actual wall clock time is even more difficult to estimate and is very 
problem dependent. For two-dimensional problems, wall clock times of 
twice the CP have been found. For three-dimensional flows, ratios of up 
to 10 could be obtained for highly irregular grids. 

The following subsections describe the input data and the output formats 
of the INTEG module. Section 5.2 describes the required inputs for the DYNDIh 
deck and Section 5.3 gives a summary or "mini- manual" of the input formats. 

A detailed discussion of each input variable is then given in Section 5.4. The 
various types of output, with examples, are described in Section 5.5. 

5.2 DYNDIM DECK 

The INTEG module of the GIM/STAR code has a "dynamic -dimens Lon" 
feature. Each dimensioned variable has a parameter in place of an actual 
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Integer. The DYNDIM program uses the UPDATE feature of CDC software 
to replace these parameters with real dimensions so that INTEG can then be 
compiled. The dimensions are under user control via input to the DYNDIM 
program. Execution of this "preprocessor" routine must be made before 
executing INTEG as explained in Section 3.4. This section describes the 
input data for the DYNDIM program. 

The required data are. input on one card with the following format: 

FORMAT (415) 


MN 


The total number of nodes in the problem. Each flow vector used in 
INTEG is then dimensioned by MN. Thus, for "small" problems, the code 
is dimensioned small and will fit in the main memory of the STAR- 100. For 
"large" problems the code is dimensioned large and the virtual memory 
handles the storage. 

IDIM 

The spatial dimension of the problem. 

IDIM = 1 axi symmetric flow 

= 2 two-dimensional planar flow 
= 3 three-dimensional flow 

This parameter is used to set the number of equations and vectors required 
for solution. For two-dimensional problems, a dimension of 1 is set for all 
variables in the third dimension. This results in a storage savings since 
the third dimension is not to be used. 
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ISPEC 


A flag denoting whether a single gas or a two-gas is to be run. 

ISPEC = 0 single ideal gas 
= 1 two ideal gases 

This parameter also controls the dimensioning of the INTEG module. If 
ISPEC = 0,' the specie continuity equation variables are dimensioned by 
1: if ISPEC = 1, these variables are dimensioned by MN. 

NSP 


The dimension of the "special node" variables. The INTEG module 
uses vectors of length MN for most of the calculations. However, for 
boundary points some of the operations are not vectorizable. These 
special cases are treated by scalar operations. The dimension of these 
scalar variables is NSP. 

The value for NSP can be set in either of two ways; 

1. NSP is output from the GEOM module. This value is subse- 
quently input to DYNDIM. 

Z. If NSP = 0 is input to the GEOM module, and the INTEG module 
then NSP = MN is the default value. 

For larger problems in which GEOM is executed in a single runstream, 
the user should set NSP from the GEOM output. If a smaller case is to be 
run, with GEOM and INTEG in the same runstream, set NSP = 0 and use the 
default value of MN. 

If a particular problem, with fixed dimensions, is to be executed many 
times, then the compiled INTEG code can be saved. On subsequent execu- 
tions then the DYNDIM program will not be needed since the relocatable 
INTEG will be called in rather than the source deck. Compiling the module 
however, does not require much CP time so that this procedure is not 
usually required. 
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5.3 INPUT CARD SUMMARY 


This subsection presents a summary of the input cards and formats 
for the INTEG module. A description of each input parameter and its options 
is detailed in Section 5.4. After a user becomes familiar with die INTEG 
inputs, this summarized input guide can be used to quickly identify each card 
and its contents. 

Three basic formats are used to input the data to INTEG; 

ALPHANUMERIC A1 
INTEGER 15 

DECIMAL ElO.O 


Integers are thus right justified in five column increments. Decimal, or 
floating point data occupy ten columns each with, preferably, a decimal point 
punched on the card (see examples. Section 8). 


Card Type 


Parameter List/Format 


1 ICASE, (ITITLE(I), 1=1, 78) 

(12, 78A1) 

2 IDIM, METHOD, ITMAX, IPRNT, ITSAVE, ISTART, 
lOTYPE, lUNITS, ITSTRT, IVISC, IDIST, ISPEC 

(1215) 

3 NN, NNX, NDX, NNY, NDY, NNZ, NDZ, NPM 

(815) 

4 DTIME, DTFAC, INCDT 

(2E10.0, 15) 

5 REALMU, REALK, GAMSl, GAMS2, WMl, WM2, DK, RK 

(8E10.0) 

6 EMU, ELAM, ERHO, ESPEC 

(4E10.0) 
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Card Type 


Parameter List/Format 


7 


8 

9 

10 

11 

12 


NNPM(I), NCPM(I), (NNCPM(I, J), J= 1. 5), 
ANGPM(I); 1= 1, NPM 

(715, ElO.O) 

(NCT(I, J,K), PXPM(I, J,K), PYPM(I, J,K), K= 1,4); 
J= 1, NCPM(I); 1=1, NPM. 

(15, 2E10.0) 

RHOZ, PZ, ASTAR, NINC, A, B 
(3E10.0, 15, 2E10.0) 

NJ, INC. NTOT. ITAN, ITYPE 
(515) 

RI. UI, VI, WI, PI, CSI 
(6E10.0) 

Nl, IC. NT 
(315) 


5.4 DESCRIPTION OF INPUT DATA 

This section presents a description of the input parameters listed In 
Section 5.3. Each parameter is identified as to its usage in the INTEG 
module with options and standard values also given. All 12 card types are 
not necessarily input for a given problem. Some of the control parameters 
on card type 2, for example, dictate which options have been selected and 
hence which input cards are required. This information is given in the 
discussion of each parameter to be input. Each input card which is read 
by INTEG is then immediately printed out to aid the user in debugging a 
problem setup. 

5.4.1 Card and Parameter Descriptions 

CARD TYPE 1 Format (12, 78A1) 


Problem Identification Labels 


ICASE 


The numerical number assigned to the given specific problem. This 
is used only for user identification of the case being run. 

ICASE > 0 identifies the case 
< 0 terminates the run 


Any number of cases can be stacked and run in sequence with the INTEG 
module. After the last case to be processed, the user can set ICASE = -1 
to allow a normal termination of the run stream. 

ITITLE 


Any alphanumeric information that the user wants to assign to identify 
specifics of the problem. The title should start in Column 3 and can extend 
through Column 80 of the card. This title is used only for identification of 
the problem. 

CARD TYPE 2 Format (1215) 

Program Option Control Parameters 


IDIM 


Flag indicating the dimensionality of the problem being run. 

IDIM = 1 indicates that the two-dimensional 
axisymmetric option is to be used 

= 2 for two-dimensional planar flow 

= 3 for full three-dimensional flow 

The only restrictions on the use of IDIM is that the input files 17 and 20 
from the GEOM module must have been generated consistent with two- or 
three-dimensional flow. IDIM is also an input in GEOM and must be con- 
sistent with its use in INTEG. 
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METHOD 


A flag Indicating whether a one- step or a two-step time integration 
method is to be used. 

= 1 one -step method 
= 2 two-step method 

The only restriction to use of METHOD is again a consistency with the GEOM 
module. If METHOD = 2, then the INTEG module expects to find two sets 
of matrices on the input file 17 from GEOM, One step methods probably 
should not be used for problems involving shock waves, or for viscous 
flows. 

ITMAX 


The maximum number of iterations to run on this execution of INTEG* 
There is no method currently in the module to stop a run other than on the 
specified number of iterations. A suggestion is to run, say, 100 iterations 
and check the results, then restart and run several hundred more steps. 

Nominal Values; ITMAX = 50 

= 500 

= 1000, etc. 

Discussions on "when a case has converged" are given in Section 7. 

IPRNT 


The print control flag. Output of the entire flow field, at selected 

th 

nodes, is made by the INTEG module every (IPRNT) iteration. For 
example. 


IPRNT =100 

provides a flowfield printout every 100 time steps. 
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A suggestion is to printout often on a debugging run and then on subse- 
quent production runs only printout every several hundred steps. 

ITSAVE 


The flowfield save flag. Output of the flow field is written on File 22 
at every (ITSAVE) iteration. For example, 

ITSAVE = 500 

will write every 500th time step out on file 22. This file can then be saved 
on permanent files or rolled out to tape. This file 22 can then be used in 
two ways: 

1. Plotting of the flow field by Module 3 is done from file 22. 

2. Restarting an INTEG run requires that file 22 be saved from 
the previous run. 

A good idea is to always save file 22 on any long run and then delete it if 
no longer needed. 

IS TART 


The restart control flag. This controls the input data sequence for 
the INTEG module. 

ISTART = 0 indicates a ”cold“ start 

run with all initial flow- 
field data read from cards. 

= N>0 indicates a restart run. 


If N > 0, then file 22 must be available for INTEG to read. The numerical 
value of N is the block number of the record on file 22 to be used as the 
starting data. 


Example: A previous run has saved file 22 at every 100 
iterations for a 1000 iteration run. A restart 
case is wanted starting at iteration 800. 
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Then N= 9 since block 9 will contain the 800th iteration. Iteration 0 is the 
first block on file 22 if the user has flagged IOTYPE<0. If lOTYPE > 0, 
then N =: 8 in the example because iteration 0 has not been written on the 
file 22. 

lOTYPE 

Indicator of the type of flowfield output to be printed. An example of 
each available type of output is shown in Section 5.5. There are three basic 
types; 

• One-line output prints only the basic flowfield variables such 
as velocity, density, pressure, etc. 

• Two-line output prints extra flow variables which are calculated 
after the calculation such as Mach number, flow angle, etc. 

• Simplified output consisting primarily of supersonic flow quan- 
tities of interest such as Mach number, pressure, etc. 

lOTYPE = _+ 1 gives one-line output 
= + 2 gives two -line output 
= + 3 gives the special purpose 

If lOTYPE > 0 
< 0 

The user should refer to Section 5.5 for examples of the various types of 
output than IN TEG can print. 

lUNlTS 

A flag indicating the units to be used on all flow variables. 

lUNITS = 1 English units are assumed 

= 2 CGS metric units are assumed 


output 


then the zero iteration 
is not printed 

the zero^ or starting itera 
tion is printed. 


5-11 



Either system may be used but all inputs must be in that consistent set. If 
dimensionless, variables are used, then set lUNITS = 2. The assumed units 
are: 


Variable 

English 

Metric 

Velocity 

ft/ sec 

cm/ sec 

Density 

Ibm/ft^ 

gm/cm 

Pressure 

Ibf/ft^ 

dyne/ cm^ 

Viscosity 

Ibm/ft-sec 

gm/cm-sec 

Thermal 

Conductivity 

Ibm/ft-sec^-R 

gm/ft-sec^-K 

Gas Constant 

ft^/sec^-R 

cm /sec -K 

Binary Diffusion 
Coefficient 

ft^/sec 

2 / 

cm /sec 

Time 

sec 

sec 


ITSTRT 


Number of the iteration at which the run is to start. This is used only 
for printout purposes to help the user keep track of the iterations on a se- 
quence of multi-iteration restart cases. For the example restart problem 
discussed in the ISTART card, the value of ITSTRT is 800. To keep all 
iteration numbers in order, the values of ITSTRT and ISTART should be 
consistent, i.e.. 


ISTART = 9 

ITSTRT = 800 

is consistent for the example problem. Again, this is for accounting purposes 
only, and an error in ITSTRT will not affect the run or the answers. 

IVISC 


Flag indicating whether the run is to be made with the inviscid Euler 
equations or the full Navier-Stokes equations. 
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IVISC = 0 Indicates an Invlscid 
case to be run 

= 1 signals a full viscous 
Navler -Stokes run 

If a viscous run Is called for, then the appropriate coefficients must be Input 
on subsequent cards. For Invlscid runs, these coefficients can be read In as 
zero since they are not used. The only restriction Is that IVlSC = 1 must be 
used If a two-gas case Is being run. Mixing of the two gas streams could 
not occur, of course, for an Invlscid flow. 

There are two viscous options available, both flagged by IVISC = 1. 

1. Real laminar viscosity can be Input as non-zero and INTEG 
will use this value. 

2. Only NDC type viscosity can be called upon. This Is used for 
numerical stability only and should not affect the problem 
answers themselves. 

Care should be taken when attempting to compute an Invlscid supersonic flow 
with shock waves. This type of problem should be run as viscous with only 
NDC type artificial viscosity In order to capture the shocks. Methods of In- 
putting real viscosity and NDC coefficients will be given subsequently. 


ID 1ST 


This Is a special Input flag to be used only for nozzle/boundary layer 

flows. 


IDIST = 0 provides an invlscid type of 
starting solution for nozzle 
flows 

= 1 provides for a viscous boundary 
layer type of starting solution 
for a nozzle flow. 


Details for Inputting a nozzle flow starting solution are given in the descrip- 
tion of card type 9. 
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ISPEC 


Flag indicating whether the problem to be run consists of one ideal 
gas or a two-gas system. 

ISPEC = 0 indicates a single ideal gas 
= 1 indicates a two-gas system 

For single ideal gases, only one continuity equation is solved. For two-gas 
systems, a global and a single species continuity equation are solved. Prop- 
erties of each of these two ideal gases are input on subsequent cards. For 
ISPEC = 0, the input of these properties can be zero, since they are not used. 
For ISPEC = 1, the values of molecular weight and specific heat ratios must 
be input. The binary diffusion for a species into the mixture is also required. 
The IVISC = 1 option must also be selected. 

CARD TYPE 3 Format (815) 

Node Point Description 


NN 

The total number of nodal points in the problem. This input must be 
consistent with the input to the GEOM module. It is read in on a card in 
INTEG rather than passed across on a file for convenience. 

I There is, in theory, no limit on NN except for the external storage 
devices of the STAR- 100 system. In practice, the number of nodes, NN 
should be selected based on accuracy and resolution of the problem. Keep- 
ing NN as small as possible, however, reduces the CP time and the large 
page faulting on the STAR system. 

NNX 


The number of nodal points in the x-coordinate direction. This param- 
eter is input (and NNY, NNZ) because of the possible irregular grid system 
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of the GIM code. To compute a time step, the INTEG module must know the 
dimensions of the grid in each coordinate direction. 

For some irregular grid configurations, die number of nodes in each 
coordinate direction may be different among the regions of the full grid. 

For this reason, we will sometimes have 

NN ^ (NNX) (NNY) (NNZ). 

As a guideline for irregular grids, the value of NNX (and NNY, NNZ) should 
be chosen for that part of the grid having the most nodes. These inputs are 
for computing a time step only. If a constant time step is input by the user, 
the values of NNX (NNY, NNZ) are not meaningful, but should not be input 
as zero . 

NDX 


This nodal parameter (and NDY, NDZ) is input to the INTEG module 
to indicate the coordinate directions in which the nodes are numbered. The 
GEOM module has the option, specified by the user, to number the nodal 
points in any order, i.e., starting along any of the three axes first. For 
INTEG to set up the vectorization scheme consistent with data provided by 
GEOM, the program must know the numbering algorithm. 

NDX is the nodal point decrement in the x-coordinate direction, i.e., 
the difference between two node numbers adjacent to one another along the 
x-axis. Consider the following two-dimensional example: 
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For tills example the input would be 

NNX = ? NNY = 4 

NDX = 4 NDY = 1 

and similarly for a three-dimensional problem. 


NNY 


Number of nodal points in the y-coordinate direction (see discussion 
of NNX) 

NDY 


Nodal point number decrement in the y-coordinate direction (see 
discussion of NDX). 


NNZ 


Number of nodal points in the z-coordinate direction (see discussion 
of NNX). For two-dimensional problems set NNZ = 1. 

NDZ 


Nodal point number decrement in the z-coordinate direction (see dis- 
cussion of NDX). For two-dimensional problems set NDZ = 0. 

NPM 


This flag controls the treatment of sharp expansion corners in two 
dimensional supersonic flow. 

NPM = 0 indicates that no special treat- 
ment of corner points is to be 
made 

= N indicates that there are N sharp 
corners to be treated explicitly 
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The input data for each corner is given on card types 7 and 8 discussed 
later. If NPM = 0, then cards 7 and 8 are not input at all. The current 
limit on the number of special corners is 10, i.e., NPM < 10 is the limit. 

CARD TYPE 4 Format (2E10.0, IS) 

Time Step Control Parameters 

DTIME 

The value (in seconds) of the time step At to be used ^ the constant 
step size option is chosen. The user has the option of using this value for 
At or having the code compute At based on the CFL condition. The next 
parameter to be input is DTFAC, which determines the option chosen. 


DTFAC 


A factor which has two functions. The first function of DTFAC is its 
multiplication by the minimum CFL number to produce the minimum At for 
each step: i.e.. 


At . 


mm 


IdTFAcI * CFL 


min 


The second function depends on the algebraic sign of DTFAC. 


If DTFAC > 0 the code will use the user DTIME 

value as long as it Is below the 

At . . 
min 

If DTFAC < 0 the code will always use the Atmin 
itself and ignore the user input 
value of DTIME. 


If the value of DTIME is ever detected to be larger than At . , then DTIME 

min 

is set equal to At . . 

^ mm 
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The CFL condition is calculated in INTEG by the following formula; 


CEL = 


lul . Ivl . Iwl 
-r— + "V" '" + ~T — + a 

Ax Ay Az 


T Ax Ay Az 


-1 


for each node point, where 


Ax, Ay, Az are the grid sizes respectively in 
the X, y, z direction 

u, V, w are the velocity components re- 
spectively in the x, y, z directions 

a is the sonic velocity. 


No computation is made for a diffusion controlled time step. The DTFAC 
parameter can be used to adjust the CFL. 


INCDT 


Iteration counter indicating the frequency at which to update the At 
calculation. The computation of all the CFL conditions for all nodes, and 
then searching for the minimum, is a costly operation. It is not usually 
necessary to do this at every time iteration. The value of INCDT can be 
used to control this calculation. For example, 

INCDT = 50 

tells the code to run with a constant time step for 50 iterations and then 
check the stability value, update the step size if necessary, and then run 
another 50 iterations. If INCDT = 1, the code will compute the CFL values 
at every iteration. If the time step is cut to satisfy the CFL, a message 
is printed to this effect. 

CARD TYPE 5 Format (8E10.0) 


Gas Thermodynamic and Transport Properties 



REALMU 


The value of the dynamic viscosity of the ideal gas system. In the 
current code, this is a constant, laminar-type value. There are no models 
in this version of GIM to compute the viscosity coefficient. The units of 
this value depend on lUNITS input on Card 2 (see this discussion). 

REALK 


The value of the thermal conductivity of the ideal gas system. This 
value is also currently a constant in the appropriate set of units. 


GAMSl 

The ratio of specific heats for ideal gas number 1. If a single gas 
system is being run (ISPEC = 0), GAMSl is the value of y for this gas. 

GAMS 2 

The ratio of specific heats for ideal gas number 2. If a single gas 
system is being run, set GAMS2=0.0. 


WMl 

The molecular weight of gas number 1. If a single gas system is being 
run, WMl is the molecular weight of this gas. The code has the option, for 
a single gas, of reading-in directly the gas constant RK or using WMl and 
the universal gas constant (see discussion of RK to follow). 


WM2 


The molecular weight of gas number 2. The code uses the universal 
gas constant and molecular weights to compute the individual gas constants 
for use in the ideal gas equation of state. 
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Ill 


DK 

The value of the binary diffusion coefficient for species 1 into the 
mixture. 


If DK ^ 0 the code will use this value in the species 
continuity and energy equations. 

DK < 0 the code will compute this coefficient as 
DK = P/l/0.7 

from the dynamic viscosity /Lt with an 
assumed Lewis number of 1.0 and Prandtl 
number of 0.7. 

For single gas systems, the value of DK is not used and can be input as 0.0. 
RK 

The gas constant for a single ideal gas. If RK > 0, the code will use 
this value for the individual gas constant. If RK= 0, the code will compute 
RK from the universal gas constant and the molecular weight WMl. 

CARD TYPE 6 Format (4E10.0) 

NDC Coefficients (See Fig. 3-2.) 

EMU (Nominal Value 0. 0 to 0.5) 

Coefficient for multiplying the computed NDC artificial viscosity. 
This, and the other NDC coefficients are used to control the amount of 
numerical diffusion in the solution. The total viscosity is computed as 
the sum of REALMU and the artificial term: 

MU = REALMU + EMU*(NDC Value) 
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ELAM (Nominal Value = -2/3) 


The coefficient to use in computing the "second" viscosity coefficient 
X in the Navier -Stokes equation is as follows; 

X ~ ELAM*4 

The value of -2/3 gives the Stokes relation. 

ERHO (Nominal Value 0.0 to 0.5) 

Coefficient to multiply the computed NDC artificial diffusion term for 
the global continuity equation. This value is usually set to zero for shock- 
free flows, but should be about 0.5 for problems which have strong shocks 
(see Fig. 3-2.) 

ESPEC (Nominal Value 0.0 to 0.5) 

Coefficient to multiply the computed NDC artificial diffusion term for 
the species continuity equation. Same usage as ERHO for global conserva- 
tion of mass. 


CARD TYPE 7 Format (715^ ElO.O) 

Node Point Information for Special 
Treatment of Sharp Corners 


Notes 


1. Input card type 7 and 8 only if NPM > 0 on Card 3. 

2. For NPM > 0, the input sequence is one card type 7 followed 
by cards type 8; then repeat the sequence NPM times. Card 7 — 
then 8, Card 7 — • then 8, etc. 

3. The number of cards type 8 to be input after each type 7 is a 
variable specified on card type 7 itself. Thus, for different 
NPM values, there may be more than one card of type 8 (see 
description of Card 8). 

4. In the subsequent input description, the subscript I runs from 
1 to NPM. 
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NNPM(I) 


Node number of the corner point to be specially treated. See the 
example below. 



Node 701 is the corner node itself. SetNNPM(I) = 701. 

NCPM(I) 

The total number of nodes, connected to the corner node, which the 
user wants to be computed using the downstream side expanded flow at node 
NNPM. In the above example, if NCPM(I) = 3, then three of the connecting 
nodes are assumed to be downstream of the Mach line from the corner. 

NNCPMd, J), J= 1, 5 

These are the node numbers of those connecting nodes which are to be 
recomputed. In the above example, NCPM(I)= 3, thus, there are three values 
of NNCPM to be input. (The J = 1,5 is for generality, simply set the re- 
maining two values to zero.) For the example, the input would be 

NNCPM(I, 1) = 801 
NNCPM(I, 2) = 802 
NNCPM(I, 3) = 702 
NNCPM(I, 4) = 0 

NNCPM(I, 5) = 0 

This then connects the corner node 701 to surrounding nodes 801, 802 and 
702. These nodes will be computed using the downstream flow from the 
corner. Nodes 601 and 602 are also connected to the corner, but are "up- 
stream" of the expansion and will be computed from the upstream unex- 
panded flow. 
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ANGPM(I) 


The expansion angle of the corner, measured from the horizontal. 



The input angle is in degrees. 


CARD TYPE 8 Format (15, 2E10.0) 

Finite Difference Information for Corner Treatment 


Notes 


1. The number of cards type 8 is 

4*NCPM(I) 

where NCPM(I) was input on Card 7 for 1=1,2,..., NPM. 

2. In the subsequent input description, the variable I runs 
from 1 to NPM, and the variable J runs from 1 to NCPM(I). 

NCT(1, J,K), K= 1, 4 

The node numbers (4 of them) which are connected to node number J, 
for treatment of corner node I. 


Note that the format is (15, 2E10.0) for Card 8. This means that the 
following card sequence is required: 


Card 8-1 NCT(I, J, 1), PXPM(I, J, 1), PYPM(I, J. 1) 

Card 8-2 NCT(I, J, 2), PXPM(I, J, 2), PYPM(I, J. 2) 

Card 8-3 NCT(I, J, 3), PXPM(I, J. 3). PYPM(I, J, 3) 

Card 8-4 NCT(I, J. 4), PXPM(I, J, 4), PYPM(I, J, 4) 
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PXPMd, J,K), K= 1,4 


The value of the "der Lvatlve-taker" matrix B for computing the 
x-derivative of node J due to the contribution of node K. These values 
are taken from printout of the GEOM module. The best method of under- 
standing this procedure is to study the example problems of Section 8. 

PYPM(1, J,K), K = 1, 4 

The value of the "derivative- taker” matrix C for computing the 
y-derivative of node J due to the contribution of node K. See discussion 
of NCT and PXPM. 

CARD TYPE 9 Format (3E10.0, 15, ZEIO.O) 

Nozzle Flow Parameter Identification 


Notes 

1. Card type 9 is input for all cases except a restart run. 

2. The data on Card 9 are actually used only if a nozzle-flow 
starting conditions type is flagged on Card 10. 

3. The data on Card 9 can be blank if it is not to be used; 
but the blank card must be there. 


RHOZ 


Value of the stagnation density of the gas flow. A one -dimensional 
isentropic equation is solved to give an initial density distribution in the 
nozzle. The flow conditions are simply expanded isentropically based on 
area ratio of the nozzle. 


P/Po = 
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PZ 


Value of tlic stagnation pressure of the gas flow. Same usage as 
RHOZ, where 

p/p. . 

AS TAR 


The area of the nozzle throat. The Isentroplc expansion, is done in 
terms of area ratio, i.e., local area to ASTAR. The code uses a subsonic 
expansion for the converging portion of the nozzle and switches to a super- 
sonic equation after it encounters an area ratio of 1.0. The test that is 
made is 


ASTAR 

A 


1 < l.OE-4 


The user should take care to ensure that ASTAR is accurate enough to allow 
the code to detect the throat of the nozzle. If this test falls, a subsonic ex- 
pansion will also be used in the diverging section of the nozzle. 


NENC 


The nodal point number increment between nodes on the nozzle center- 
line and the upper wall. This is used to obtain the coordinates of the nodes 
to compute the local area. In general, NINC is the number of nodes in the 
y (or radial) direction minus 1. If there are 20 nodes radially in the nozzle, 
then NINC =19. 

A, B 


These are parameters that the user can set to obtain a "boundary- 
layer -like" velocity distribution in the nozzle. These parameters are used 
if IDIST = 1 was input on card type 2. The equation used is 

k = A+ B*x (xis the axial coordinate) 
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then 


M = M (1 “ y/y 

o 'w 

where 

M = local Mach mxmber at location (x, y) 

= local isentropic Mach number at the centerline 
of the nozzle (x, y= 0) 

y = local radial coordinate of the node 

X = local axial coordinate of the node 

y^ = radial coordinate of the upper wall of the nozzle 

This gives a value of M^ at the centerline, M = 0 (no slip) on the nozzle wall, 
and a "boundary-layer-like" distribution in between the boundaries. 


CARD TYPE 10 Format (515) 

Nodal Designators for Initial Conditions Input 


Notes 


1. Any number of cards type 10 may be used to initialize a flow 
field. All nodes can be input on a single card or each node 
can be input on a separate card. The usual case is some- 
where between these two extremes. 

2. If ITYPE = 0 is input on card type 10, then a card type 11 must 
immediately follow. If ITYPE X 0, no card 11 is input. 

3. A -1 card (Columns 4-5) must be input as the last card in a 
type 10 sequence to terminate rea^^ng of initial conditions 
data. 

4. Card type 10 (and 11 if required) are not input for a restart 
case unless changes are being made (see Section 5.472)” TRe 
-1 card must be present even on a restart case. 




Node number of the first nodal point to be initialized by this card 
type 10. (NJ = -1 terminates the input of Card 10.) 
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INC 


Node number increment to NJ to be used for inputting a sequence of 
nodes on one card. Set INC = 0 if only one node is to be set by this card. 
If INC > 0, this will be added to NJ to allow the user the convenience of 
incremental inputs. 

NTOT 


The total number of nodes to be set by this card type 10. 

Example; NJ = 600, INC = 2, NTOT = 5 will cause nodes 600, 602, 604 
606 and 608 to be initialized to the same values (input on Card 1 1 to follow). 


ITAN 


A flag which can be used to allow the user to input either; (1) velocity 
components, or (2) total velocity. If 

ITAN = 0 the code assumes that the user 

will input u, V, w velocity compo- 
nents . 

= 1 the code will expect a total velocity 

to be input. The individual compo- 
nents will be computed by the code 
using the flow angle information 
from the Geometry file. 


For boundary conditions which require inviscid free slip to be maintained, 
ITAN = 1 should be input to ensure that accuracy is maintained in applying 
the "tangency" condition at walls. This option can also be called upon to 
start the flow in the proper direction to enhance the convergence to steady 
state. 


ITYPE 


The critical flag that tells the code which type of initialization is to 
be done. 
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ITYPE 


= 0 allows the user to input the initial 
conditions on cards (see card 
type 11). 

= 1 signals the code to use the nozzle 

flow initialization option (see card 
type 9). 

= 2 sets the initial conditions for this 

card type 10 input to the same values 
as the previous card type 10 node 
values. 

= 3 allows the user to code a subroutine 
USERIP in order to initialize the flow 
field. This subroutine is provided in 
skeleton form with the code. Often a 
flow problem can be best initialized 
from the output of another code, from 
interpolated data or another source. 
This option is provided to give the user 
total flexibility in initializing the flow 
field. The reason being that the initial 
data determines to a great extent, the 
amount of computer time it takes to 
converge to steady state. (Examples 
are given in Section 8 of the use of 
USERIP inputs.) 


CARD TYPE 11 Format (6E10.0) 

Flowfield Initial Conditions 


Notes 

1. Card type 1 1 is input following a card type 10 only if 
ITYPE = 0 is input on card type 10. 

2. If ITYPE = 0 on Card 10, then a type 11 card must follow. 

3. Cards type 1 1 is not input on a restart case unless changes 
are being made (see Section 5.4.2). 

RI 


Mass density to be assigned to the nodes specified on the correspond- 
ing card type 10. 
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UI 


Velocity component in the x-coordinate direction to be assigned as 
initial conditions. If ITAN = 1 was input on the corresponding card type 10, 
the UI is to be the total velocity rather than the x-direction component. 

yi 

Velocity component in the y-coordinate direction. If ITAN t= 1, then 
set VI = 0.0 and the code will compute the component from geometric inputs. 

WI 

Velocity component in the z-coordinate direction. If ITAN = 1, then 
set WI = 0.0 as is the case with VI. For two-dimensional problems, input 
WI = 0.0. 


PI 


Static pressure to be assigned to the nodes designated on the corre- 
sponding card type 10. 


CSI 


Mass fraction of gas species number 1 to be assigned as initial condi- 
tions. This must be between 0.0 and 1.0. For a single gas case, IS PEC = 0, 
input CSI = 0.0. 


CARD TYPE 12 


Nodal Output Control 


Notes 


1. The nodal point flow field can be output in any order for con- 
venience of the user. The nodes do not have to be printed in 
sequential order. 


5-29 



2. The maximum number of nodes to be printed in NN, the total 
number in the problem. A given node can be printed more than 
once, if it is convenient, but the total number of output stations 
must not be more than NN. 

3. Any number of cards of type 12 may be input. For example, 
each node to be printed can be specified on a separate card, or 
all nodes to be printed can be specified on a single card. The 
usual case is somewhere between those two extremes. 

4. A -1 card (Columns 4-5) must be last card in the sequence of 
type 12 cards. This causes the code to start the execution 
cycle. 


m 

The node number of the first nodal point on this card to be printed. 
For example, N1 = 698, then node 698 will begin the printing controlled by 
this card. 

IC 


The node number increment to N1 at which printing is desired. For 
example, N1 = 100, IC = 3 will cause output of nodes 100, 103, 106, 109* . . . 


NT 


The total number of nodes to be printed by this card type 12. For 
example, N = 200, IC = 2, NT = 5 will cause output to be printed for nodes 
200, 202, 204, 206 and 208. 

Example of Card Type 12 Sequence: 


1 

1 

100 

101 

2 

50 

501 

5 

200 

-1 
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5.4.2 Special Input Sequence 

There are always "special" cases to be Input with most any flowfield 
code. There are often problems which can best be initialized in some non- 
standard way. There is no known method of single card input which will 
cover all possibilities. In the GIM/STAR code, we have provided two options 
for handling these cases. The RSTART option allows a problem to be initial- 
ized from the GIM solution of another case. This is described in the subse- 
quent discussion. The second choice is to use the USERIP option. This 
allows the user to write a special purpose subroutine to initialize a problem. 
At Lockheed we have found this to be the most efficient way of initializing 
many problems. The better the initial guess of the flow field, the faster is 
the convergence to steady state. The USERIP option is also discussed in 
the following pages. 


RSTART Option 


The RSTART feature of the code has two basic functions; 

1. It allows a problem to be "restarted" from a previous run. If 
a case does not converge for some reason, the problem can be 
restarted from some point with different parameters. This also 
allows a case to be run in segments. A run can be made for, 
say 200 iterations and then stopped. The run can then be re- 
started at iteration 200 and run again. 

2. It allows parametric studies to be made efficiently. A case for 
a given set of flow conditions, (pressure, Mach number, etc.), 
can be converged for a configuration. The output can then be 
used to initialize a case for different pressure, Mach number, 
etc. This permits steady state results to be obtained in much 
less computer time than a "cold start" run. 

To use the RSTART option, the user proceeds as follows; 

1. On the first run of a sequence, set ITSAVE > 0 and save file 22 
on a tape or permanent file. 

2. On the second run of a sequence, set ISTART > 0, and copy the 
previously save file onto file 22. 

3. Repeat this sequence for each restart case. There is no limit 
to the number of times a case can be restarted. 
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An important item to remember on a restart run is that the initial 
conditions cards type 10 and 11 are not input unless changes are to be made. 
For the standard RSTART usage, only the -1 card of sequence types 10 and 11 
is needed. The option does exist for allowing the user to restart a run and 
change some of the conditions if desired. To use this option, the following 
sequence should be used: 

1. Set ISTART > 0 and copy the saved file onto file 22. 

2. Input cards type 10 and 11 in the same manner as a ’’cold start” case. 

3. Conclude the type 10 and 11 sequence with the standard -1 card. 


The example problems given in Section 8 show an example of the use of the 
RSTART option. 


USERIP Option 


This option allows the user to write a special purpose subroutine to 
initialize a flow field. The skeleton of this subroutine is provided with the 
GIM/STAR code. To exercise this option, the following procedure is used: 

1. On card type 10, set the parameter ITYPE = 3. The subroutine 
USERIP is then called with arguments, NJ, INC, NTOT, IRTN. 

The first three of these are just the inputs from the card type 
10. (See this section for definitions.) The parameter IRTN 
should be set by the user in USERIP. If IRTN= 0, the code pro- 
ceeds to store the flow variables in the standard input manner. 

If IRTN = 1, the code will assume that the user has stored all 
data and return is made to read another card type 10. 

2. The second step is to code the required data into USERIP. This 
is done using the UPDATE feature of the CDC software. The 
ID ENT is as follows: 

*IDENT USER 

*D PAGE2.142 

User Codel 


The variables that must be set in this subroutine are RI, UI, VI, WI 
and PI as defined in the discussion of card type 11, If IRTN= 0 is set, then 
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the above variables must be returned to the INPUT subroutine. If IRTN= 1 
is set, the code expects the user to set the following variables; 

RHO(I), UVEL.(I), VVEL(I), 

WVEL(I), P(I), ENER(I) 

The example problems given in Section 8 show an example of the 
USERIP option. To use this particular option of the code efficiently requires 
some experience and practice with its use. An experienced user will, how- 
ever, find this to be a very efficient means of inputting the important initial 
conditions. 

5.5 OUTPUT DESCRIPTION 

The output of the INTEG module consists of two types; 

• Printed flowfield data, and 

• Output file on which the flowfield data 
are stored. 

The output file (22) is a formatted file used in the GIMPLT module to plot 
contours of the flow field (see Section 6) and for restarting a run. The 
printed output from INTEG consists of the following types; 

• The input c.ards are printed out to allow the user to check the 
data. 

• At each time iteration, a printout is given of the unsteady 
derivative parameters. 

• At user selected iterations, the flow field itself is printed. 

Each of these types of output is now briefly described. The reader should 
refer to the corresponding figures while reading these descriptions. 

5.5.1 Printout of Input Cards (Figs. 5-2 and 5-3) 

The first page of output from the INTEG module is a summary of the 
input cards. The variable names and nomenclature for this printout are the 
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Fig. 5-2 - Example Printout of Input Cards 
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Fig. 5-3 - Problem Summary Output 


same as described in the INPUT guide (Section 5.4). The user should verify 
that the intended input got into the code correctly. This page is followed by 
problem summary statements which give the code's interpretation of the 
options requested by the user. Again, this should be checked to ensure that 
the proper options were selected. The final line on the input summary is 
titled: 

"Unsteady Derivative Reference Values." 

The quantities are defined at current iteration "zero" and are the sum of 
squares of the unsteady derivatives; 

NN 

E (^i)^ 

i=l 

The order of this printout is 

p^, (pu)^, (p’v)^, (pw)^, (PE)^, (pC)^ . 

. 2 

For two-dimensional problems, the (pw) is not printed and for single gas 

. 2 . 

cases, the (pC) is not printed. The values of these variables are different 

for each problem and set of initial conditions, but generally are of the order 
5 20 

of 10 to 10 . The primary use of this data is for debugging of a new problem 

set-up and to check for convergence of a set of restarted runs (see Section 7). 

5.5.2 Sum of Squares of Unsteady Derivatives (Fig. 5-4) 

The output philosophy of the GIM code is to print "some information at 
all iterations and all information at some iterations." At each time step, the 
current values of the "sum of squares" of the unsteady derivatives is printed. 
This information has proved to be very useful in determining convergence of 
a case. The printed values start at 1.0 since they are normalized by the 
reference values at iteration zero; In general, these derivatives should start 
to go down as the solution gets closer to steady state. In some cases, they 
may go up for a while as the code is adjusting the solution. Convergence is 
generally obtained after the values drop several orders of magnitude and then 
stop changing. 
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The iteration number, ITER, is followed by the current value of the 
minimum time step, DTMIN. The normalized sum of squares of the unsteady 
derivatives are then given for each variable in the U vector. A header is 
printed at each new page to identify the derivative being printed. See Section 
7 for further comments on the unsteady derivative method for converging a 
case. 


5.5.3 Flowfield Output (Figs. 5-5 through 5-9) 

Printout of the flow field is given at user specified iterations. In addi- 
tion, the user can also select the node points to be printed. These are con- 
trolled by flagging the appropriate input cards. A large volume of paper can 
be generated if all nodes are printed for a three-dimensional problem. Care- 
fully selected nodes could be printed during the convergence of a case, and 
then all nodes printed after steady state is reached. 

The format of this output depends on two input parameters: 

IDIM; lOTYPE 


If IDIM= 2, only those variables pertinent for two-dimensional flows are 
printed. If IDIM= 3, the full three-dimensional flow field is output. The 
value of lOTYPE determines the amount of printout at each node point. 


IQTYPE = 1 : This will be called the "one line" output and consists of 
the basic calculated flow parameters; density, velocity, energy and pressure. 
Figures 5-5 and 5-6 are examples of this output for two- and three-dimensional 
flow, respectively. The nomenclature for the printout is as follows: 


NODE 
X, Y, Z 
RHO 
U, V, W 
QVEL 
E 


— node point number 

— Cartesian coordinates of the NODE 

— density 

— velocity components in x, y, z directions 

— total velocity (u^ + v^ + w^) 

— total energy 
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Fig. 5-5 - Two-Dimensional One Line Output 
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Fig. 5-7 - Two-Dimensional Two Line Output 
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11CRAI IHN NHMHEK 1 10 

DTIHEx 
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05 TlMt 
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, 20 no() 060 t +02 

.0 
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1 76 
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.6576751 5E+03 

.21 120307E+01 

.200001>6uc. *02 

.0 


4 

l'*2 

,2752446be-02 
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,21 1 288966*01 

.20000060* *02 

.0 
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?06 
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.0 
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.0 
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Fig. 5-9 - Special Output Format 
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p 

IB 


— pressure 

— node point type (see Section 4) 

lOTYPE = Z; This will be called the "two -line" output and consists of 
the basic flow variables plus auxiliary calculations. Figures 5-7 and 5-8 
are examples of a two-line output for two- and three-dimensional flow, re- 
spectively. In addition to the basic variables printed out for lOTYPE = 1, 
the following are given for lOTYPE = 2. 


GAM 

~ local value of the ratio of specific heats 

CS 

— mass fraction of gas species number one 

PH IX 

— flow angle in the x-y plane 

PHIY 

— flow angle in the y-z plane 

SOS 

— local sonic velocity 

T 

— temperature 

M 

— Mach number 

If a two-gas code is being run, then lOTYPE = 2 will give a printout of the 
species continuity equation solution. 

lOTYPE = 3: 

This will be called the "simplified output" option and 

consists of only auxiliary flow variables; 

RHO 

— density 

P 

— pressure 

M 

— Mach number 

PH IX 

— flow angle in x-y plane 

PHIY 

— flow angle in y-z plane 

IB 

— node point type 


This option is generally useful for simple cases such as inviscid supersonic 
flow where only pressure and Mach number are wanted. It can also be used 
during the debugging stage to reduce the volume of numbers to be viewed. 
Figure 5-9 is an example of this type of output. 
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At the top of each new page, the current iteration number and time 
step are also printed. Following each flowfield iteration is a STAR run time 
printout. The CP time in seconds and the wall clock time in seconds is given. 
The values are the elapsed times since the last printout. For exarnple, if 
output is requested every 10 iterations, and the CP time = 1.62 then each 
iteration is taking 0.162 sec on STAR. 

The units for the output variables are the same as the input set of units. 
This is set by the user input value of lUNITS (see Section 5.4). 


5.5.4 Error Messages 


The current INTEG module has a few error diagnostic messages that 
are printed if the code detects an irregularity. This is not a complete set 
of possible errors, but have been found to be most common ones which occur. 
A brief description of these seven error messages follows. 

1. TOO MANY OUTPUT NODES REQUESTED, MAXIMUM = NN 

This message indicates that the user has made an 
error on the node printout card type 12. The maxi- 
mum number of nodes that can be printed is NN, the 
total number of nodes in a problem. 


2. TIME ITERATION DID NOT CONVERGE 

This message is printed if the solution "blows -up" 
before ITMAX iterations are reached (see Section 7). 


3. UNSTEADY DERIVATIVES GREATER THAN TOLERANCE 

The normalized sum of squares of the unsteady 
derivatives are monitored by the code. If any of 
the values goes up by more than 5 orders of magni- 
tude, this message is printed. It usually means that 
the solution is about to "blow-up." The time step 
parameters, DTFAC, INCDT are a possible source 
of the trouble. The NDC coefficients should be 
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checked and possibly adjusted. An error in the initial 
conditions can often cause this nnessage. 


4. INITIAL CONDITIONS INPUT FOR NODE NUMBER GREATER 
THAN NN 

A common error is the input of initial conditions data for 
a non-existent node number. If this error message is 
printed, the user should check the input cards type 10. 


5. INVALID lOTYPE 

The value of lOTYPE, on the input card type 2, must be 
either 1, 2 or 3. 


6. INVALID BOUNDARY CONDITION TYPE 

The INTEG module checks the boundary condition flags 
IB which are input from the GEOM module. All nodes 
must have IB values between 0 and 9. The user should 
check the GEOM module input data deck. 


7. TWO-GAS CASE REQUESTED WITHOUT VISCOUS OPTION 

If a two-gas problem is being run, the INTEG module 
requires that the IVISC = 1 option is selected. This is 
done to ensure that the mixing coefficients for the equa- 
tions are all consistent. The GIM formulation does not 
use an "inviscid slip-line" approach for shear flows and 
hence must use the viscous terms to mix a flow from two 
different streams. If ISPEC = 1, then the user must also 
set IVISC = 1. 


Other possible error messages may be printed by the STAR system 
itself. The user should refer to their documentation for explanation of 
these diagnostics. 
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6. GRAPH ICS MODULE GUIDE (GIMPLT) 


6. 1 USE OF THE MODULE 

The GIMPLT module may be used to generate plots of the finite element 
grid, velocity vectors, and/or pressure, temperature, density and Mach number 
contours. 

Portions of the flow field may be plotted separately to obtain enlarged 
views of certain regions, such as corner regions, throat sections, etc. These 
regions may be plotted to the same scale, if desired, so they can be placed 
together to form a large plot of the entire flow field. 

The description of the region to be plotted including its node/element 
topology, view angle, etc., is defined by a "Plot Specification." Once these 
"Plot Specifications" are input, grid plots, velocity vector plots, and/or con- 
tour plots may be generated for any, or all, of the flowfield regions described 
by the set of "Plot Specifications." 

Many of the input data parameters to the GIMPLT module have default 
values which are used if the input parameter field is left blank or set to zero. 
Most of these default values will be satisfactory in most cases, relieving the 
user of some data input. 

The GIMPLT module may be executed following the IN TEG module for 
generation of velocity vector, contour, and/or grid plots. If grid plots only 
are desired, the GIMPLT module may be executed following the GEOM module. 
When checking out a new model, it is recommended that plots of the grid be 
generated and examined before proceeding with INTEG executions. 
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The file requirements for the GIMPLiT module are illustrated in 
Fig. 6-1. 



Fig. 6-1 - File Requirements for GIMPL.T Module 


The GIMPLiT module resides on the STAR Access Station (CDC 6400 
"Z" machine). Hence, all files needed must be available on the Access 
Station when executing the GIMPLT module. Since the GEOM and INTEG 
modules execute on STAR, files 20 and 22 must be transferred back to the 
Access Station via the TOAS command. See Section 3 of this report and/or 
the STAR Programming Manual for details. 

Following the execution of the GIMPLT module, the system PLOT 
post-processor must be executed to obtain the actual plots. Several 
graphics devices are available at Langley. Refer to the Graphic Output 
System User's Guide for details. 

6.2 INPUT CARD SUMMARY 

This section presents a summary of the input cards and formats for 
the GIMPLT module. A detailed description of each input parameter, its 
options, and default values (where present) is given in the following 
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Section 6,3. After a user becomes familiar with the GIMPLT input param- 
eters this input card summary can be used to quickly identify each card and 
its contents. 

Three basic formats are used to input the data to the GIMPLT module; 

ALPHANUMERIC Axx 
INTEGER IS 

DECIMAL ElO.O 


Note that some fields are skipped for some card types as indicated by nX in 
the format description. 


Card Type 


Parameter List/Format 


1 ITITLE(l), ITITLE(2) 

(2A40) 

2 NX, ITERAD, ITRBLK, KDIM, ISP 

( 515 ) 

3 GAMMA, FACTOR, RK, PO, TO, RHOO 

(6E10.0) 

Specs . 

S-1 NPLT, STITLE, IVIEW, ISYM, ITHETl, lAXISl, 

ITHET2, IAXIS2, IXTABL, lYTABL, VFAC 

(15, 5X, A20, 815, ElO.O) 

S-2 NTYPE, JO, IJUMP, JJUMP, NI, NJ, IPRNT 

(715) 

Grid 

G-1 'GRID', lOPT, ICSCLE, NSPECS, (ISPEC(I), 1= 1, 

NSPECS) 

(A4, IX, 15, 25X, 215, 715) 

G-2 (ISPEC(I), 1=8. NSPECS) if NSPECS > 7) 

(45X, 715) 

VVEC 

'VVEC, lOPT, NITER, ICSCLE, NSPECS, (ISPEC(I), 
1=1, NSPECS) 

(A4, IX, 215, 20X, 215, 715) 
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Card Type Parameter LLst/Format 

VVEC 

V-2 (ISPEC(I), 1=8, NSPECS) (if NSPECS > 7) 

(45X, 715) 

I-l (ITER (I), 1=1, NITER) 

(1615) 

Contours 

C-1 ITYPE, lOPT, NITER, NC, ITABLE, INCR, ICSCLE, 

NSPECS, IS PEC (I), 1=1, NSPECS) 

(A4, IX, 515, 5X, 215, 715) 

C-2 (ISPEC(I), 1= I, NSPECS) (if NSPECS > 7) 

(45X, 715) 

I-l (ITER (I), I = 1, NITER) 

(16I5) 

L-1 (CVAL(I), 1= 1, NC) 

(8E10.0) 


6.3 DESCRIPTION OF INPUT DATA 

This section presents a description of the input parameters listed in 
Section 6.2. Each parameter is identified as to its usage in the GIMPLT 
module with options of each shown. 


Example plots are given for a few selected cases. All possible plot 
combinations cannot be illustrated due to large numbers of possibilities. 
The user can generate those options of interest for the example problems 
as a learning tool for understanding the GIMPLT module. Finally, certain 
restrictions which apply to GIMPLT on the STAR- 100 system are given. 

CARD TYPE 1 Format (2A40) 

Problem Identification 
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ITITLEd), ITITLE(2) 


This problem identii ication title appears on contour and velocity vector 
plot frames as two lines of forty characters each. This title identification 
should be made consistent with the problem identification for the correspond- 
ing INTEG run which generated the flow field to be plotted. 

CARD TYPE 2 Format (515) 

Problem Size and Type Description 

NX 

The total number of nodal points in the model which is being plotted. 

NX must be consistent with the GEOM and INTEG runs which generated the 
grid and flow field. 

ITERAD 

The iteration number of the first set of flowfield information on the 
INTEG File 22 which is being used as input to the GIMPLT module. This 
parameter may be left blank if GRID plots only are being generated. 

ITRBLK 

The number of iterations per block on File 22. This corresponds to 
the iteration save increment ITSAVE used in the INTEG run which created 
File 22. 

This parameter may be left blank if GRID plots only are being gen- 
erated. 

KDIM 


The dimensionality of the problem. 

KDIM = 2 for 2-D planar or axisymmetric 
KDIM =3 for 3-D 
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ISP 


The two -gas flag used in INTEG. 

ISP = 0 for single gas 
= 1 for two gases 

This parameter is needed to define File 22 content and may be left 
blank if GRID plots only are being generated. 

CARO TYPE 3 Format (6E10.0) 

Gas Constants and Reference Values 
(This card may be left blank if contour plots are not being generated.) 
GAMMA 

The ratio of specific heats as input in the INTEG module. 

FACTOR 


This is a unit conversion factor. 

Set FACTOR = 32. 174 if lUNITS = 1 in INTEG 
Set FACTOR =1.0 if lUNITS = 2 in INTEG. 


RK 


The gas constant as input in INTEG. 

PO, TO, RHOO 

These are reference values which are used to normalize pressure, 
temperature and density, respectively, for contour plots. These may be 
stagnation values if desired. If contour plots of the unnormalized quantities 
are desired, PO, TO and RHOO may be set to 1.0. 
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A set of "Plot Specifications" follows. Each "Plot Specification" 
defines the node/elcment content of the plot, its orientation, view angle, etc. 

A description of the four-node elements which constitute the flowfield region 
represented by each plot specification is required. These element outlines 
always appear on GRID plots but normally do not appear on velocity vector 
or contour plots unless requested by the user. 

A set of two-node line co.nnectors is available to describe the outline 
of the flowfield region and/or any distinguishing features, such as centerlines, 
flow boundaries, etc. Once this library of "Plot Specification" is established, 
GRID plots, velocity vector, or contour plots may be generated for any or all 
of the flowfield region by selecting the proper "Plot Specification." The input 
data requirements for a single "Plot Specification" follow. The data for sev- 
eral "Plot Specifications" may be stacked. A card containing a - 1 in Columns 
4 and 5 terminates the "Plot Specif Lcation" data library. 

CARD TYPE S-l Format (15, 5X, A20, 815, ElO.O) 

Plot Specification Labeling and View Angle Options 

NPLT No Default 

Plot specification identification number. NPLT may be any positive 
number less than or equal to 100. The plot specifications need not be iden- 
tified in ascending order. NPLT = -1 terminates the plot specification data. 


STITLE No Default 

The plot specification title. This title will appear on grid plots Iden- 
tifyliig the plotted region, e.g., CORNER REGION. 

IVEEW Default IVIEW = 3 

The view axis which is normal to the plotted plane. 
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ISYM 


Default ISYM=0 


A symmetric reflection parameter. If ISYM^O, the plane normal to 
axis ISYM is a symmetry plane. The synunetric reflection about this plane 
is plotted. This parameter is normally left blank. 

The next four parameters on this card, ITHETl, lAXISl, ITHET2, 
IAXIS2 along with the view axis, IVIEW, are used to select a non-standard 
axis orientation on an oblique view angle. 

ITHETl Default ITHETl = 0 

Rotation in degrees about axis lAXISl. 

lAXISl No Default 

Axis about which first rotation is performed. 

ITHET2 Default ITHET2=0 

Rotation in degrees about axis IAXIS2. 


IAXIS2 No Default 


Axis about which second rotation is performed. The standard view 
angle consists of the z-axis (3) directed toward the viewer (out of the paper) 
with the X and y axes directed to the right and top respectively, as illus- 
trated below. 








+z 


X 


IVIEW = 3 
(Standard) 


+x 




IVIEW = 1 


+y' ^ 

IVIEW = 2 


z 


If a different axis orientation is desired, or if an oblique view angle is 
required, rotations may be performed by using the ITHETl, lAXISl, ITHET2, 
and IAXIS2 parameters. Following are two examples. 
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-Z j ►x 

y 

IVIEW = 3 
ITHETl = 180 
LAXISl = 1 

Note that the second rotation (ITHET2, IAXIS2) is performed about the ro- 
tated reference frame resulting from the first rotation. 

When creating contour plots, an option exists for labeling the contours 
and constructing a label ID table. (See discussion of contour plots later in 
this section.) This table is normally placed in the upper left-hand corner of 
the plot frame. The following parameters may be used to position the ID 
table elsewhere on the plotted frame, 

IXTABL Default IXTABD=20 

X raster value for upper left-hand corner of label ID table. 

lYTABL Default IYTABL=980 

Y raster value for upper left-hand corner of label ID table. 

1023, 1023 
Plotted Frame 


0, 0| 1 1023,0 




X 

IVIEW = 1 

ITHETl = 30 1AXIS1 = 2 
ITHET2=-30 IAXIS2=3 
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VFAC 


Default VFAC =1.0 


A length factor used to lengthen or shorten the vectors on velocity 
vector plots. (See discussion of velocity vector plots later in this section.) 

Examples 

VFAC= 2.0 doubles the length 

VFAC = 0,5 halves the length 

CARD TYPE S-2 Format (715) 

Plot Specification Node/Element Content 

Card type S-2 is used to specify the node/element content of a plot 
specification. Any number of cards of this type may be input to describe 
a single plot specification. The last card for a particular plot specification 
should be blank except for a - 1 in Columns 9 and 10, which terminates the 
data for that plot specification. 

NTYPE Default NTYPE=0 

Element generation control parameter, 

NTYPE = 0 a single element is input 

- Z a network of two -node line 
connectors is generated 

= 4 a network of 4-node elements 
is generated 

JO No Default 

For NTYPE = 2 or 4, JO is the starting node in the generated network. 
For NTYPE = 0, JO is the first connected node. JO= -1 terminates the se- 
quence of card type S-2 input for a plot specification. 
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IjUMP 


No Default 


th 

For NTYPE = 2 or 4, IJUMP is the nodal increment in the i direction 
for the element network generation. For NTYPE = 0, this is the second con- 
nected node, i.e., node JO connects to node IJUMP. 

JJUMP No Default 

fVfc 

For NTYPE = Z or 4, JJUMP is the nodal increment in the j direction. 
For NTYPE = 0, JJUMP is the third connecting node, i.e., node JJUMP is con- 
nected to node IJUMP. 

NI Default NI = 1 

For NTYPE = 2 or 4, NI is the number of elements generated in the i^ 

direction. For NTYPE = 0, NI is the fourth connecting node, i.e., node NI is 
connected to nodes JJUMP and JO. 

NJ Default N J = 1 

th. 

For NTYPE = 4, NJ is the number of elements generated in the j di- 
rection. If NTYPE = 2, NJ is the number of lines of connecting nodes in the 
fVi 

j direction. For NTYPE = 0, NJ is not used. 

IPRNT Default IPRNT = 0 

IPRNT =0 no printout 

= 1 nodal connections are printed 
for each element generated by 
this card 

For NTYPE = 2, a network of 2-node line connectors are generated as illus- 
trated below. 

Where the 'origin* and 'terminus' of a line are connected to nodes J1 
and J2, respectively. The program generated NI*NJ lines, as follows. For 
1=1, through NI, and J = 1 through NJ. 

Jl = JO + (J-1)*JJUMP+ (I-1)*IJUMP 
J2 = Jl + IJUMP 
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Example; NT YPE = 2, JO = 6, I JUMP =12, JJUMP = 2, NI = 5 , NJ = 5 


J = 1 
J = 2 
J = 3 
J = 4 
J = 5 



For NTYPE = 4, a network of 4~node elements are generated as illus- 
trated below. 

Where Jl, J2, J3, and J4 are the nodes connected by the element, 
the program generates, fox 1=1 through NI, and J = 1 through NJ, 

Jl = JO + (I-l) «1JUMP + (J-1) * JJUMP 

J2 = Jl + IJUMP 

J3 = J2+ JJUMP 

J4 = J3 - IJUMP = J1+ JJUMP 

Example: NTYPE =4, JO = 6, IJUMP = 2, JJUMP = 12, NI = 4, NJ = 5 



J=1 J=Z J=3 J=4 J=5 
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The remaining data describes the content and form of; 


finite element grid plots, 

— velocity vector plots, and 
~ P, T, p and Mach number contour plots. 

Any number of finite element grid (GRID), velocity vector (VVEC), and/or 
contour (PRES, TEMP, DENS, MACH) definition decks may be Included in 
any order. 

A card containing STOP in Columns 1 through 4 (or a blank card) 
terminates the entire GIMPLiT input data sequence. 

CARD TYPE G-1 Format (A4, IX, 15, 25X, 215, 715) 

Grid Definition Card 


grid No Default 

The first entry on this card consists of the alphameric characters 
GRID in Columns 1 through 4. This field defines this card as a grid defi- 
nition card, calling for grid plots to be generated. 


lOPT Default lOPT = 0 


Option parameter for optionally numbering the nodes and/or shading 
elements on grid plots. 


lOPT 


0 no numbering or shading 

1 nodes will be numbered 
Z elements will be shaded 

3 nodes numbered and elements shaded. 
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ICSCLE 


Default ICSCLE = 0 


Constant scale parameter which permits plot specifications to be plotted 
to the same scale such that portions of a flow field appearing on separate plot 
specifications may be plotted on separate frames and placed together to form 
the complete flowfield grid. 

ICSCLE = 0 all plot specifications are plotted to 
the maximum size on each frame 

= 1 all plot specifications called for on 
this GRID card will be plotted to the 
same scale. 

NSPECS Default NSPECS = 0 

The number of specifications to be plotted as a result of this GRID 
card. If NSPECS is input as zero, or left blank, all plot specifications will 
be plotted. If NSPECS > 0, NSPECS values will be read on the remainder 
of this card and continuing on following card(s) if necessary. 

ISPEC(I), 1=1, NSPECS) 

List of plot specifications to be plotted as a result of this card. If 
NSPECS > 7, continue on subsequent card(s) starting in Columns 46 through 
50. 


CARD TYPE C-2 Format (45X, 715) 

List of plot specifications is continued on cards of this type if 
NSPECS > 7. 

CARO TYPE V-1 Forniat (A4, IX, 215, 20X, 215, 715) 

Veloeity Vector Card 
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VVEC 


No Default 


The first entry on this card consists of the alphameric characters 
VVEC in Columns 1 through 4. This field defines this card as a velocity 
vector parameter card, calling for velocity vector plots to be generated. 

lOPT Default lOPT = 0 

Option parameter for optionally superimposing the grid on velocity 
vector plots. 

lOPT = 0 grid will not be superimposed 
= 1 grid will be superimposed. 

NITER Default NITER = 1 

Number of iterations for which velocity vector plots are to be generated. 

ICSCLE Default ICSCLE = 0 

Constant scale parameter. Refer to GRID card for description. 

NSPECS Default NSPECS = 0 

The number of specifications to be plotted as a result of this VVEC 
card. If NSPECS is input as zero, or left blank, all plot specifications will 
be plotted. If NSPECS > 0, NSPECS values will be read on the remainder 
of this card and continuing on following card(s) if necessary. 

(ISPEC(I), 1=1, NSPECS) 

List of plot specifications to be plotted as a result of this card. If 
NSPECS > 7, continue on subsequent card(s) starting in Columns 46 through 50. 

CARD TYPE V“2 Format (45X, 715) 

List of plot specifications is continued on cards of this type if 
NSPECS > 7. 
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CARD TYPE I”! Format (16I5) 


(ITER (I), 1=1, NITER) 

List of iteration numbers for which velocity vector plots are desired. 

If more than l6 iterations are desired, continue on subsequent cards of this 
type. If plots of initial conditions only are desired, NITER may be set to 
zero, or left blank, and Card 1-1 omitted. Iteration 0 will be plotted auto- 
matically. 

CARD TYPE C-1 Format (A4, IX, 515, 5X, 215, 715) 

Contour Plot Card 

ITYPE No Default 

The first entry on this card identifies the type of contour plots desired. 
Four types of contour plots are permitted, identified by the following four- 
character alphameric names. 


PRES - 
TEMP - 
DENS - 
MACH - 

Pressure contours (P/PO) 
Temperature contours (T/TO) 
Density contours (RHO/RHOO) 
Mach number contours 

Default 

lOPT = 0 


Option parameter for superimposing the grid on contour plots. 

lOPT = 0 grid will not be superimposed 
= 1 grid will be superimposed. 


NITER Default NITER = 1 

Number of iterations for which contour plots are to be generated as a 
result of this card. 
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NC 


Default NC = -10 


Number of contour values to be plotted as a result of this card. NC 
may be input as a positive or negative number. If NC is positive, the actual 
contour values are read in on card types L-l. If NC is negative, the con- 
tour levels are calculated automatically based on the maximum and minimum 
values for the flowfield variable for which contours are being plotted. The 
absolute value of NC must be less than or equal to 50. 

IT^ABLE Default ITABLE=0 

This is a contour labeling flag. If ITABLE = 0, the contours are not 
labeled. If ITABLE/ 0, the contours are labeled and a table is constructed 
showing the actual contour values. 

If ITAB1 jE= 1, the table will be placed directly on the contour plot 
frame at raster locations IX, lY specified on card S-1. 

If ITABLE = -1, the table will be placed on a separate frame at raster 
locations IX and lY. This is sometimes necessary if there is not enough 
room on the contour frame itself for the table. 

INCH Default INCR = 7 

This is a contour labeling increment controlling the labeling density 
on the contours. It is approximately the number of grid points between 
contour labels. 

If the contours are labeled too often, increase INCR. If the contours 
are not labeled often enough, decrease INCR. 

1CSCL.E Default ICSCLE = 0 

Constant scale parameter. Refer to GRID card for description. 
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NS PECS 


Default NSPECS = 0 


1 


The number of speclflcatLons to be plotted for this set of contours. 

If NSPECS is input as zero, or left blank, all plot specifications will be 
plotted. If NSPECS > 0, NSPECS values will be read on the remainder of 
this card and continuing on following card(s) if necessary, 

(ISPEC(I), 1=1, NSPECS) 

List of plot specifications to be plotted for this set of contours. If 
NSPECS > 7, continue on subsequent card(s) starting in Columns 46 through 50. 

CARD TYPE C-2 Format (45X, 715) 

List of plot specifications is continued on cards of this type if 
NSPECS > 7. 

CARD TYPE I“1 Format (l6l5) 

(ITER (I), 1=1, NITER) 

List of iteration numbers for which contour plots of this type are de- 
sired. If more than 16 iterations are desired, continue on subsequent cards 
of this type. If contour plots of Initial conditions only are desired NITER 
may be set to zero or left blank, and Card I- 1 omitted. Contour plots for 
iteration 0 will be generated automatically. 

CARD TYPE L-1 Format (8E10.0) 

(CVAL(I), 1= 1, NC) 

List of values for which contours will be plotted. IF NC was input as 
a negative number, omit this card. If NC is greater than 8, continue on 
subsequent cards of this type. 
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6.4 OUTPUT DESCRIPTION 


Printed output from the GIMPLiT module consists essentially of just 
an echo printout of the input parameters. When input parameters are left 
blank the default values will be printed out. Additional printout can be called 
for on card type S-E by setting IPRNT = 1. This causes the nodal connections 
to be printed out for each element generated by that S-E card. Otherwise, 
just the generation parameters as input will be printed. In addition, follow- 
ing the echo print of each GRID, velocity vector, and contour plot control 
parameters, a message is printed indicating the number of plot frames 
generated for that command. 

The main output from the GIMPLT module are the plots. Figure 6-E- 
is an example of a GRID plot for a lEl-node duct case. This figure shows 
the two types of plots, with and without nodal numbering. Figure 6-3 is ex- 
amples of solution contours for the lEl-node source flow expansion. Shown 
are velocity vectors and Mach contours at steady state (approximately EOO 
iterations). 

Figure 6-4 is an example of a GRID plot for a converging-diverging 
nozzle. The model was only for half of the symmetric nozzle, but the 
GIMPLT module can plot the reflected half as shown. Figure 6-5 is a Mach 
contour plot for this nozzle problem showing the details of flow that can be 
resolved. Figure 6-6 is an example of a three-dimensional grid generated 
for an expanding duct. Many combinations of plot types can be obtained. 

These figures illustrate the type of capability available. 

6.5 RESTRICTIONS 

There is a limit to the size problem which can be plotted with the 
GIMPLT module in its present form. The GIMPLT module resides entirely 
on the 6400 "Z" machine which does not have the paging capability of the STAR. 
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Fig. 6-4 - Grid Plot for Nozzle Flow 
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Fig. 6-5 - Example Mach Contours for InviscLd Nozzle Flow 
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Fig. 6-6 - Three-Dimensional Grid for an Expanding Dact 
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The total core requirements for plotting including the program storage, sys- 
tem library, and all data arrays are given by the following formula: 

KFLjq = KMAXjq + 16000jq 

where 

KMAXjq = 5040 + 8#NX (NX = number of nodes) 

KMAX is the dimension of the array *'A'' in the GIMPLT MAIN program which 
is normally set to 15000 and must be changed to plot a larger model. KMAX 
is also a program variable which must be set in the MAIN program. 


For KMAX = 15000, the maximum number of nodes in a model to be 
plotted is 


NX 


15000 - 5040 
8 


1245 


The field length requirement for KMAX = 15000 is 

KFLjo = ISOOOjo + l 6000 jo 
= 31000jo 

= 74430g 

The CM parameter on the JOB card must be at least as large as KFLg. A 
field length of ZOOOOOg (65536^^) would permit the following. 

KMAXjq = 65536jg - IOOOO^q = 49536 

NX = 49536 - 5040 _ nodes. 

O 
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7. USING THE GIM/STAR CODE 

The GIM/STAR code described herein is operational on the STAR- 100 
system at NASA -Langley. A number of example cases have been success- 
fully computed by Lockheed -Huntsville through a remote batch terminal. As 
with any large fluid dynamics code^ the GIM-ON-STAR version may experience 
problems for some flow cases. It is recommended that the user contact the 
orginators at Lockheed -Huntsville to discuss any such problems which may 
arise. The following subsections summarize a number of items that should 
prove useful in running the code on the STAR system. 

7.1 SETUP PROCEDURE 

A first step, that we have found useful before setting up a problena on 
GIM, is to have an open discussion of the case with people knowledgeable in 
that type of flow. Major problem areas can.be isolated, the type of flow 
patterns expected can be identified and modeling techniques discussed. The 
next step is to set up a geometry model of the flow domain, decide on the 
number of grid points needed and run the GEOM module. Execution of the 
GIMPLT module should be made next, if the grid is relatively complex. A 
plot of the grid, especially for three dimensions. Is often times the only way 
to get it right. We usually require more than one try to get the nodal distri- 
bution into the code properly. 

The next step is to set up the INTEG module input deck. Selection of 
appropriate initial conditions is critical for converging a case in a reasonable 
amount of run time. It is usually a good idea to plot contours of the initial 
conditions and to check iteration zero before making a long run. Always 
make a short run (10 iterations or so) and check the unsteady derivative 
printout, time step and spot check a few critical nodes. We have found that 
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the most efficient way to run a case is to use the "RSTART" option. Run 
a few hundred iterations and save File 22, Check the flow field, either by 
printout or plots and restart the case running a few hundred more. As the 
problem converges to steady state the unsteady derivatives should go down 
several orders of magnitude and then stop changing. 

The final step is to execute the GIMPLT module to give the desired 
flow field contour maps. This sometimes requires several runs to get the 
exact pictures that you want. An additional item that we ususally do after 
a case has converged, is to run INTEG for a "dummy restart" case and 
simply printout the entire flow field of nodes. This will provide a permanent 
record for the files of the complete problem. 

The following discussion describes some of the procedures that we have 
found useful in making the code work. 

7.2 GRID SELECTION 

The choice of mesh sizes for the problems that can be solved with GIM 
is dictated primarily by the truncation error of the difference equations and, 
in some cases, by the stability criteria. The only way in which a final mesh 
size may be chosen with complete confidence is to run the problem with suc- 
cessively smaller mesh sizes until little or no change is observed in the 
results. Such a procedure may be impractical from the standpoint of com- 
puter time required, and in this case an alternative, although somewhat less 
desirable, solution is to run the problem for three mesh sizes and then extra- 
polate to zero mesh size. This method will not produce satisfactory results 
if the mesh sizes chosen are not fairly close to the size required to give a 
reasonable solution. 

No quantitative statements can be made about the choice of mesh sizes 
for general domains. Smaller mesh sizes are usually required in regions 
of more rapidly changing velocity and pressure. This means that a fine 
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«niesh is required close to the wall in boundary layer problems, as well as 
I close to leading edges. In the region close to a wall, or in the mixing region 
I of streams of different velocities, it is usually desirable to employ a finer 
' mesh size than that needed far out in the free stream where flow gradients 
are less steep. 

In general, the use of different mesh sizes in the same flow field 
does not significantly change the GIM method of solution. This is true since 
a special form of the difference equations is used at the point of mesh size 
change. GIM uses ’'general interpolation" functions to account for mesh size 
changes. Whenever possible, the use of a grid with rapidly changing mesh 
size should be avoided. It is better to stretch the mesh and have uniformily 
changing mesh size rather than abrupt changes. 

7.3 INITIAL CONDITIONS 

If only a steady state solution is to be computed, the initial conditions, 
theoretically, do not effect the results. In practice there are three primary 
ways that initial conditions effect a GIM code run. 

1. Convergence to steady state can be done with much less 
computer time if the initial conditions are chosen as close 
as possible to the expected solution. 

2. Improper initial conditions, such as violating the ideal gas 
law, mass conservation, etc., can cause spurious behavior 
of the relaxation solution. 

3. Sharpdisconti ny i ties in the initial conditions, such as a 
"sharp" shocl^an cause early instabilities to set up and 
ultimately diverge the solution. 

Thus the user should take special precautions to start a problem as 
best he can. In the GIM/STAR code, a special input option called USERIP 
is provided which will allow a user to set up initial conditions for special 
circumstances. We recommend using this option and did so in both the 
two-dimensional and three-dimensional Scramjet example cases. Another 
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alternative Is to start with the sharp discontinuities and use a very sniall 
time step until the discontinuity gets washed out and then switch to a realistic 
step size. 

7.4 DIFFERENCE OPERATORS 

One of the advantages of the GIM formulation is the flexibility of the 
finite difference algorithm. We have experimented somewhat with various 
algorithms and have found that for explicit calculations, the two-step 
MacCormack scheme works best for a wide range of flows. If the user has 
any doubt about the choice of scheme, we recommend this MacCormack-type 
algorithm In this fully elliptic, explicit version of the GIM code. 

For certain cases, a choice of two steps in the MacCormack algorithm 
must be made. For most applications to flows with shock waves, we have 
found that "forward-forward-(forward)" and "backward-backward-(backwardy 
operators work best. For supersonic, inviscid flows without shocks, we have 
had some success with "backward-forward-(forward)” for step 1 and ’’backward- 
backward-(backward)’' for step 2. This scheme always differences '’backward" 
into the supersonic flow and does not allow the downstream side to cause any 
influence (as physically it does not). 

A simple one-sided, one-step, Euler scheme has also worked well for 
certain subsonic flows. The truncation error is of first order and generally 
is unacceptable for flows with any appreciable gradients. It is suggested that 
the beginning user print out selected nodal analogs in the GEOM module and 
check them to see that the desired scheme has indeed been input properly. 

7.5 CONVERGING A CASE 

A definite well formulated and programmable criteria for determining 
when a Navier-Stokes solution has converged is not known to these authors. 
Consequently, in the GIM/STAR code we have left some of this burden on the 
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user. We have, however, provided some assistance in making these decisions. 
The code prints out the " sum of squares of the unsteady derivatives" at each 
iteration as explained in Section 5.5. The user should monitor these quantities. 
Theoretically, they should go to zero as actual steady state is reached. Due 
to finite discretization of a problem, this will of course not be true. As steady 
state is approached, the "sums square" should go down. Based on our exper- 
ience, they should go down several orders of magnitude and then stop changing 
significantly. If the solution is going unstable, these "sums square" will begin 
to increase in value. The code is currently set to stop the calculation if any 
of the values go up by five orders of magnitude. 

Another technique that we commonly use is to plot contours at selected 
iterations and overlay the plots on a light table. Significant changes can be 
detected and the solution Iterated for more steps. If the contours exhibit no 
change for approximately 100 steps or so, we generally stop the calculation. 

A careful look at the full printout of the flow field can then be made to detect 
ahy special area of flow which may not have settled out. Any other technique 
known by the user can, of course, be applied to the GIM solution as with any 
finite difference code. 

7.6 ESTIMATING RUN TIME 

A definite formula for run times of the GIM/STAR code has not been 
established at this time. The GDC 7600 version has such a formula, but 
due to the vectorization on STAR, this formula no longer applies. The follow- 
ing tables gives the CP times for the GEOM and INTEG modules for the two 
example problems which were computed on STAR. 




GEOM CP 

INTEG CP 

Config. 

Nodes 

(sec) 

(sec for 100 steps) 

2-D Scramjet 

940 

27 

10 

3-D Scramjet 

7904 

1009 

124 
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Several notes should be made about this table. For the two-dimensional 
Scramjet problem, the INTEG module converged in ~ 1000 steps which re- 
quired about 100 seconds of CP time, compared to 27 sec total for the GEOM 
module. If a number of parametric cases were to be run, then the GEOM 
module time would be insignificant compared to the INTEG times. For 
larger three-dimensional cases, such as the 7904-node Scramjet, the run 
time would be essentially comparable among the two modules, for one case. 

Still a 1000 sec job on STAR is not insignificant as with the smaller problems. 

The actual CRU time of a run is even more difficult to estimate. For 
problems which require very few large page faults, the CRU ’’time” of a 
GIM/STAR run is essentially all CP time. However, for problems such as 
the 7094-node Scramjet cases, the large page faulting on STAR does cause 
the CRU time to be significantly higher than the CP. At this writing, no 
definitive criteria have been established to relate the run time (cost) with the 
CP seconds. It depends largely, perhaps nonlinearly, on the way the grid is 
constructed, and the manner in which the LOAD card is arranged. Our recom- 
mendation at this time is to make a short INTEG run on the problem and deter- 
mine the number of large page faults. Extrapolation can then perhaps yield a 
run time estimate for the particular case. 

7.7 TROUBLESHOOTING 

If a GIM/STAR code run '’bombs,” there are a number of items which 
should be checked early in the debugging stages. Among these are: 

• Check the dynamic dimension inputs to ensure that the 
code has been properly compiled. 

• Check the grid plot for proper connectivity and distribu- 
tion of points. 

• Analyze the nodal analogs to ensure that the finite differ- 
ence scheme is a "good" one. 

• Hand calculate a time step at a selected set of nodes to 
ensure that violation of the CFL is not being inadvertently 
done in the INTEG module. 
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• Increase the NOC coefficients if the INTEG solution 
appears to be slowly blowing-up. They can be reduced 
again after the solution has settled down. Nominal 
values of the NDC coefficients are *^0.1 to 0.5. We 
have used values as high as 1.5 to stabilize a trouble- 
some case. 

• Double check the initial conditions for violation of the 
physics of the problem and for existence of unusually 
sharp gradients. 

• Check that the boundary conditions, flags and values, 
are correct as the user intended. 


During the initial usage of GIM/STAR, version SE-1, problems may occur 
which are not covered in this manual. If this happens, please call the authors* 
to discuss a particular set of difficulties. 


*(205) 837-1800. 
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8. EXAMPLE PROBLEMS 


Five example data setups are presented in this section to Illustrate 
the input formats described in this report. Full details of each case are 
not given since the purpose here is to show the user the proper way to input 
the cards. The input guides (Sections 4 , 5 and 6) should be consulted while 
reading this section. 

The five cases consist of: 

1. Supersonic Source Flow (2-D) 

2. Nozzle Flow (2-D) 

3. Scramjet Shear Layer Flow (2-D) 

4. Scramjet Shear Layer Flow (3-D) 

5. Internal Duct Flow (3-D) 

Each of these cases is now described briefly followed by a listing of the input 
cards for each module of GIM/STAR. 

8.1 SUPERSONIC SOURCE FLOW 

This simple two> dimensional supersonic source flow case is used to 
test the code on new installations. It may prove useful for the beginning 
user of GIM to follow this simple example as a first case. The configuration 
is shown in Fig. 8- la and consists of a 15 deg expansion in a planar duct. 
Source-like circular arcs are used as grids to render the flow effectively 
one -dimensional. Figure 8- lb shows the input data cards for DYNMAT, 
GEOM and DYNDIM modules. Figure 8-lc gives the input cards for INTEG 
and GIMPLT modules. 
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**•* SUPERSONIC SOUrICt flow two 0IM£NS10NAU PLANAR •»*■*» 

»* DYNMAT JoiPOT ** 

121 2 


** GEOM Input ** 
2-0 souwce FLOW CASE 


1 

2 2 

0 

0 


o 

0 1 




1 *0 
] 

l.OOOOOE 

-a 

l.OOOE- 

16 l.OOOoE-8 l.OOOE 

1 

2 1 

2 



4 

6 4 

0 


1 

1 1 

1 1 0 

□ 

0 


0*0 

0.0 


0.0 

0.0 

0*0 

0.0 


0.0 

0.0 

o.s 

0.0 


0.0 

0.0 0*0 

I .0 

0.0 


0.0 

C.O 0.0 

0.96593 

0.25882 


0.0 

15*0 0.0 

O.Afi296 

0.1294 1 


0.0 

15.0 0.0 

1 1 

1 0 

0 



1 

1 121 





1*0 


-I 


121 


** DYNoiM Input ** 
2 0 0 


1 •OOOOOE-S 


Fig. 8-lb - Supersonic Source Flow Example Case (DYNMAT, 
GEOM. DYNDIM) 
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*•** SUPeWSOMIC SOUr^CE FLOW TWO OI«»l£NbI ONAL PLANAM 



** INTtiG 

input ♦* 


1 

2-0 SOuPCf 

FLOW CASE 


2 

2 200 

2'JO 200 

0 -3 

J21 

11 11 

1 1 1 

0 0 

O.SOOOE-2 1.0 

1 


0.0 

0.0 

1.4 

O 

. 

o 

0.0 

0.0 

0.0 

o 

. 

o 

1 *o 

1.0 

1.0 

0 

1 

2 121 

1 0 


1 .4 

2.0 

o 

. 

o 

o;o 

- 1 




1 

1 121 



- 1 




-1 

END OF 1 





** 

GIMPlT 

INPUT 

** 


2- 

-O SOURCE FLOW CASE 



121 

0 

200 

2 

0 



1 .4 


1 .0 

1. 

0 


1 *0 

1 


12 I -NODE SOURCE 

FLOW 

3 

4 

I 

1 

Li 

10 

10 

0 

2 

1 

1 

llO 

10 

2 

0 

0 

1 

1 1 1 

D 

0 

0 

o 

0 

1 1 

121 

0 

0 

0 

0 


- 1 






- 1 







GRID 

0 

0 

0 




GRID 

1 

0 

0 




VVEC 

o 

2 

v) 




0 

2o0 






PRES 

0 

2 

-10 

1 



0 

2q0 






MACH 

0 

2 

-10 

1 



0 

2oO 






STOP 








2 0 0 0 0 
O 

0.0 0*0 0*0 1.0 

• O 0*0 

1*0 1*0 


l«0 1*0 

0 0 1 0 2 0 0 1*0 


Fig. 8-1 c - Supersonic Source Plow Sac^mple Case (IN^TEG, GIMPLT) 
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8.2 NOZZLE FLOW CASE 


Next in the degree of complexity for internal flow is a transonic/ 

supersonic flow in a converging -diverging nozzle. Figure 8 -2a is a sketch 

of the Space Shuttle Main Engine (SSME) nozzle contour used in the example 

calculations. The engine has been scaled down here for example only. The 

-3 

throat radius is 5 x 10 cm. The equation shown on this figure is the wall 
contour used in GEOM to generate the nozzle geometry. Figure 8 -2b is an 
input set-up for DYNMAT, GEOM and DYNDIM for this nozzle contour. Fig- 
ure 8 -2c shows the corresponding INTEG module data set-up. The GIMPLT 
module input is given in Fig. 8 -2d to plot a grid and pressure/Mach contours 
from an INTEG calculation. 

8 . 3 TWO-DIMENSIONAL SCRAM JET SIMULATION 

The flow field analyzed involved mixing the exhaust from a two- 
dimensional Sc ramjet afterbody nozzle with freestream. The problem con- 
figuration and the flow properties of the two flow streams are shown in Fig. 
8-3a. The flow streams have different values for the ratio of specific 
heats y. Details of this problem are given in Ref. 1. Figure 8-3b shows 
the input cards for the DYNMAT, GEOM, and DYNDIM modules. This case 
is initialized in the INTEG module using the USERIP subroutine explained 
‘ in Section 5. Figure 8 -3c is a FORTRAN listing of this subroutine for the 
2-D Scramjet problem. The INTEG input cards are given in Fig.8-3d for 
the USERIP optional input. The GIMPLT module input cards are given in 
Fig.8-3e. Finally, Fig. 8-3f is an example INTEG data deck for the 
RSTART optional input. It is this type of data deck that can be used to 
converge a case to steady state. 

8.4 THREE-DIMENSIONAL SCRAMJET SIMULATION 

A three-dimensional shear layer, resulting from the interaction of a 
nozzle exhaust stream with the free stream, both beneath and beside the 
nozzle, was computed using the GIM code. The configuration, shown in Fig. 
8 -4a consists of a rectangular nozzle suspended below a body or wing. 
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1.97893 - 0.77040726 ♦(Z - 0.235910) - 33.525936 



Fig. 8-2a - Nozzle Contour Used in Example Case (SSME Scaled Down) 


**** FLOW IN A Two D I i-^tNS 1 UNAu NOZZLE **** 


** OYNMAT IriPoT 
19 Z 


** GEuM ti'^PUT ** 

TWO JlMcNSIONAL Wv^ZZi-t EXAMPLE CASE 
I 2 Z 0 0 

O 0 V 

l.O i*OOE-0o 1*0E-16 l*Ot-Ob l«0t-16 l*Ot-Ob l.O l.OE-Od 

c 

111 11212 I 

4 9 A 0 0 O 1 

2T 21 O 0 21 0 

0*0 0*075 0*ib 0*<:2b 0*30 0*375 0*44 0*50 

0*3b 0*60 0*65 0*70 0*75 0*80 0*833 0*867 

0*90 0*925 0*95 0*975 1*0 

-9.84i4t-3 3*60’->0E-3 0*0 0*0 

0*0 I.OOOOE-2 0.0 0.0 


-1 .500E-2 

0.0 0.0 

O 

» 

O 



-9.04 14E-3 

0*0 0*0 

0*0 

0.0 

0*3 

-3.000E-3 

0.0 0.0 

0.0 

0*0 

0*70 

0.0 

0.0 0.0 

0.0 

0.0 


0.0 5.00000E-3 0.0 

0*0 

0*0 


-9*04 14rl-3 

8.6050E-3 0*0 

0*0 

0*0 

• 3 

-7.6950E-J 

0.121OE-J 0*0 

-25*417 

0*0 

0*42 

-2.1460E-3 

S.48J9e~3 0*0 

-25*417 

0*0 

0.70 

-1 .bOOCE-2 

8.6050E-3 0*0 

0*0 

0*0 


I 1 1 

27 1 




4 8 

4 9 0 

0 1 



66 2 1 

0 0 21 

0 



O 

. 

O 

6*960005-3 0*0 

o 

• 

o 




5, COOOE-O J.23591 102.1846 1123*9917 - • 77o*072b-o2 .44 765b 

2.80000E-3 0*0 0*0 0.0 0.0 .09 

5.0000H-2 0*0 0*0 0*C 0.0 

5.C000E-2 3.C463E-Z O.C 18.710 0.0 

1.17950E-3 5.3946E-J 0*C 37.0 0.0 .09 

21 I 0 0 

1 1 10 J 

- 1 


** DYNJIM IhPuT ** 
1932 0 0 


Fig. 8 -2b 


Nozzle Example Case (DYNMAT 


GEOM, 
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♦* INT£G INPUT ** 


2 Two OI MLfMbl w-viflL nozzle tXAi'^.PLE CASE 



2 

50 

50 

o 

o 


92 

21 

2i 

1 0 

• ■oOOO”- 

-9 

1 .0 


loo 

0.0 

.2 

• 6f'*t30L 

-4 

0.0 

-.67 

1 .01 JOE+6 

1.5354 

. 2 

y .OOOOu-3 

1 

1 

2 1 

1 

1 

2? 

1 

2 1 

1 

1 

4 J 

1 

21 

1 

1 

64 

1 

21 

1 

1 

bb 

1 

2 1 

1 

1 

] 06 

1 

2 1 

1 

1 

1^:7 

1 

2 1 

1 

1 

146 

1 

21 

i 

1 

16V 

1 

2 1 

i 

1 

I VO 

I 

c J 

1 

1 

1 

1 

21 

1 

1 

232 

1 

21 

1 

1 

2b 3 

1 

21 

1 

1 

274 

1 

21 

1 

1 

295 

1 

21 

1 

1 

31 6 

1 

21 

i 

1 

337 

1 

21 

i 

1 

3b6 

1 

21 

1 

1 

379 

1 

21 

i 

1 

400 

1 

21 

1 

1 

4ir 1 

1 

21 

i 

1 

442 

1 

21 

1 

1 

4o3 

1 

2 1 

1 

1 

404 

1 

2 1 

i 

1 

50b 

1 

21 

i 

1 

'-■•2 6 

1 

2l 

i 

1 

‘_)4 7 

1 

21 

1 

1 

bo8 

1 

21 

i 

1 

5B9 

1 

2 1 

1 

1 

63 0 

1 

2 1 

1 

1 

bJ 1 

1 

21 

1 

1 

6b2 

1 

2 I 

1 

1 

o / j 

1 

£. 1 

1 

1 

tV4 

1 

21 

1 

1 

715 

1 

2 1 

i 

1 

736 

1 

2 1 

1 

1 

7o7 

1 

21 

I 

1 

778 

1 

2 1 

i 

1 

799 

1 

21 

i 

1 

820 

1 

21 

1 

1 

84 1 

1 

21 

1 

1 

862 

1 

21 

1 

1 

8b3 

1 

21 

1 

1 

904 

1 

21 

i 

1 

925 

1 

21 

1 

1 

946 

1 

2 1 

1 

1 

9o7 

1 

2 1 

1 

1 


-2 2 0 1 0 0 

0 O 

0*0 0*0 U «0 0*0 

0.0 

20 0»0 0»0 


Fig. 8-Ec - Two-Dimensional Nozzle Example Case (INTEG) 
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1 

21 

1 

1 

1009 

1 

21 

1 

V 

lOJO 

1 

21 

1 

1 

JOS J 

1 

21 

1 

1 

1072 

1 

21 

1 

1 

1093 

1 

21 

1 

1 

1 1 lA 

1 

21 

i 

\ 

I i JS 

1 

21 

1 

1 

1 1:^6 
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21 

i 

1 

1 1 77 

i 

21 

1 

1 

} 198 

1 

21 

1 

1 

1219 

I 

21 

i 

1 

I2A0 

1 

21 

1 

1 

1261 

1 

21 

1 

1 

J2o2 

1 

21 

i 

1 

1 303 

1 

21 

X 

1 

132A 

1 

21 

i 

1 

134 5 

1 

21 

1 

1 

1366 

1 

21 

1 

1 

1367 

1 

21 

1 

1 

JA08 

1 

21 

1 

1 

1 A29 

1 

21 

i 

1 

1 ASO 

1 

21 

1 

1 

J A 71 

1 

21 

i 

J 

1 A92 

1 

21 

■ 1 

1 

1513 

1 

21 

1 

1 

J53A 

1 

21 

i 

1 

IbSS 

1 

21 

1 

1 

1S76 

1 

2 1 

i 

1 

1597 

1 

2) 

i 

1 

161G 

1 

21 

1 

J 

1639 

1 

21 

X 

1 

1660 

1 

21 

1 

1 

1 661 

1 

21 

1 

1 

1702 

1 

21 

i 

1 

1723 

1 

21 

1 

I 

1 7AA 

1 

21 

1 

1 

1765 

1 

21 

L 

1 

1 7d6 

1 

21 

1 

1 

iao7 

1 

21 

1 

1 

1828 

1 

21 

i 

1 

1049 

1 

21 

i 

1 

1870 

1 

21 

1 

1 

189 1 

J 

21 

1 

1 

1912 

J 

21 

1 

1 

~ 1 





1 

1 

1932 



- 1 





2 end 

OF 

rtDN 




Ftg. 8-2c (Concluded) 
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** GIMPLT livPUT ** 


Tuto dimensional nozzle example case 

}932 O “bO Z O 

1*535 I .0 Z.07a5b+7 l,0lo0t+6 1*0 1*0 

1 TH'^'jAT RtoIoN 


a 

1 

1 

21 20 

45 

£- 

1 

21 

20 4b 

2 


1 

21 




- 1 




2 

AFT END 

2E6I0N 


4 

946 

1 

21 

46 

2 

94 6 

21 

20 46 

2 


I9l2 

1932 




- 1 




3 


complete nozzle 


4 

1 

1 

21 20 

91 

2 

1 

21 

20 91 

2 


1 21 
1912 19J2 
~ 1 


- J 


GRID 

O 






1 

2 

1 

2 



GRID 

0 






0 

1 

3 




GRID 













PRES 

o 

1 

24 

-1 

J 7 


0 

1 

3 




50 













• Ol 15 

0 

.0125 


0.013 


0.0175 


AJ 

O 

. 

o 


O *045 

0.03 

0.035 

• 04 


• 045 


.03 


• 06 


.07 


O.Oti 

• 1 

• 2 

0-3 


0*4 


0*b 


0*6 


0.7 


o.b 

.05 

.90 

MACH 

0 

1 

26 

-1 

17 


o 

1 

3 




50 













• J 


• 4 


.6 


. 6 


1.0 


1 .2 

1 *4 

1 « u 

1 *o 


2.0 


2.3 


2.4 


^r.5 


2.6 

2.7 

^ . 6 

2*'y 


3*0 


3.1 


3.2 


3.3 


3.4 

3.5 

0.6 

3.7 


3.3 











VVEC 

0 

1 





0 

1 

3 




50 













VVEC 

0 

1 





1 

2 

1 

2 



30 













MACH 

1 

1 

0 

1 

17 


0 

1 

3 




50 













PWaS 

1 

1 

0 

1 

1 7 


o 

1 

3 





50 

STOP 


Fig. 8~2d - Two-Dimensional Nozzle Example Case (GIMPLT) 
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i-ateral Distance, y/y 



^ig.8-3a 


Two-Dimensional Shear Flow 


Configuration 
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** OYNMAT IrtPUT ** 
940 2 


** GEOM Input ** 


TWO dimensional 5 C 1 =^AMJET case 


2 

2 2 

0 1 




0 

0 1 





1 • 0 

1 . 0 £- 

08 l.OE 

-16 1 * 0 £ 

-08 1 * 0 £- 

16 l * 0£-08 1*0 

2 






I 1 1 

1 1 

11 




24 

9 4 

0 0 

0 1 



23 

16 0 

0 16 

0 



0*0 

0.05 

0*10 

0 * 15 

0*20 

0*25 0*30 

0*40 

0.45 

0*50 

0.60 

0*70 

0*80 0*90 

0*2 

0.21021 

0*0 

0*0 

0*0 


1*0 

0.21021 

0*0 

0*0 

0*0 

0*363636 

1.2 

0 . lUVlti 

0*0 

- 6.0 

0*0 

0*454545 

3 * 0 

0*0 

0*0 

- 6*0 

0*0 


3 . 0 

2 • 302 1 2 

0*0 

' 20*0 

0*0 


0*2 

1 *283 

0*0 

20*0 

0*0 


0*2 

1 *0 

0*0 

0*0 

0*0 

0*4 

1 

1 1 

1 




29 

8 4 

9 0 

0 1 



21 

16 0 

0 16 

0 



6 .C 

0.0 

0*0 

0*0 

0*0 


6 . 0 

1 

1 

3*39403 

0*0 

20*0 

0*0 


1 1 

1 




tt 

8 19 

0 0 

0 1 



£ 1 

13 0 

0 13 

0 



0.0 

0*15 

0*30 

0*4 

0*5 

0*565 0*625 

0.75 

0*8125 

0*875 

0.9375 

1 *0 


3.0 

- 2.0 

0*0 

0.0 

0*0 


6.0 

- 2.0 

0*0 

0*0 

0*0 


6.0 

0.0 

0*0 

0*0 

0*0 


3.0 

0.0 

0*0 

0*0 

0*0 


16 

1 0 

1 





I 1 100 

- 1 


1 * 0 E -06 


0*35 

1*0 


0*6875 


** OYNOIM INPUT ** 
940 2 1 3 ia 


Fig.8-3b - Two-Dimensional Sc ramjet Case (DYNMAT, GEOM, DYNDIM) 
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**«* update -^UD to Ui>E •UbEHIP* SUBROUTINE **** 
*♦ Td^o dimcnsional schamjet case ** 


* I DENT USER IP 
*D PAGE2* 14£ 

C 

C ** USER ROuTlNt TO INPUT 2-D SCRAMUET CASE ** 

C 

DIMENSION Ap( 10*4) « UMAX! 4) tAETA<10*4»«DAr(4«10)»PE(4)»PP(4) 
R£A0<5«5B0) ( UMAX ( J ) « J> 1 *4 ) « UNUM 

5b0 FORMAT (SI 5) 

DO 550 J=I *4 
JMX* JMAX( J) 

REA0(5«59( ) ( API JJ« J) « JJ«I ONX) 

Rcl AD ( 5 • ) I AtTAIUU* J) «UU=1 »uMX) 

550 continue 

590 FORMAT (BEI 0.4 ) 

DO 560 JsltjNUM 

READ (5f 595) I DA T I JU « J ) « J J= 1 « 4 > 

560 continue 

595 FORMAT <4E I 0,4 ) 

XLlNEsO.O 
PP( 1 ) =AP{ 1 . I ) 

PP( 2) =AP( 1 ) 

PP( 3) =AP( UMAX ( 4) *4) 

PP( A ) *AP( UMAX I 3 ) *3) 

DO BOO 1=1 »nTOT 
Wl =0.0 

itan=o 

IF(1,LT.6Q9) 00 TO 600 
C ** FREc. stream conditions ** 

Wl =0*003223 
Ul=5791 *0 
VI =0.0 
PI=105.B 
CSI =0.0 
GO TO 700 
600 continue 
I 1=1/16 
I I=I-I 1*16 

IFclI.EU.l) Xt-lNEsXLlNE+l .0 
IF( I l.EU.O) 1 1=16 
ETA=(XL1NE-1 .0)/42.0 
y I = 1 1 

EPS=( Yl-l.Oj/15.0 
CHI =ETA 
DO 620 J= I «4 
IFCJ.GT.2) CHI=EPS 
1F(CH1.GT. 1,0) CH1=1.0 
JMx= JMAX( J) 

DO 6 10 L.=2.JMX 
E2 = L 

IF(CHl.l_E«AETAILfU) ) GO TO 615 

610 continue 

615 continue 
C1=U2-1 

PSI =(CHI-ActAI L l.J) )/( AETA(L2»U)-AETA(L1 «U) ) 

PE( U)=AP(El , J)+PS1*( AP(L2.U)-AP<L1 * J) ) 


Fig. 8-3c - Two-Dimensional Scramjet Case (INTEG Input by USERIP) 
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620 CONTINUE 

EPS1»1*0-EPS 

Pl=< 1.0 -ETA)*(Pc( 3)“EPS1*PPC I )-EPS*PP(A) ) +ETA* ( PE ( 4 ) -EPS I *PP ( 2 ) 
S-EPS*PP(3) )+EPSl*PE( 1 )+EPS*PE(2) 

DO 650 J=2«jNJM 

1F(PI .Lt^DATlAt J) > GO TO 660 
650 continue 
660 continue 

J1 s JJ- 1 

NKR=(PI“DAT(4*J1) >/(DAT{4* JJ)-0AT<4* Jl ) ) 

NI=DAT( l« Jl )+i^Rf<*(DAT ( 1»JJ)-DAT( l«Jl I ) 

UlsDAT(2t J1 )+KRR*IDAT(2t JJ)-DAT(2* JI ) ) 

Vi=DAT(3t^l )+r<RK*(OAT( 3« JJ ) -DAT ( 3* J| > ) 

UI =SQHT<UI **2+Vl *«2) 

Cbl=l*0 

1F(XLINE»GT,1*0) ITAN=1 
IF(EPS.GE« 1 .0) ITANsl 
IF(EPS*LE#0,0) ITANal 
700 continue 

C ** CHECK FOR TANGENCY ** 

IF (1 TAN .EQ. 0) GO TO IS 
THxYxTHETA( I ) /57. 29578 
THxZ=0.0 

call TANGNTtOI t THXY«THXZ*UI «Vl «Wl ) 

15 continue 

c ** STOWE Initial conditions ** 

WHO ( 1 > =w I 

UVEL( I ) =Ul 

VVEL( I ) =vi 

P( I ) =PI ^factor 

OVELsSOWT < UI **2 + VI**2 > 

IF(1SPEC .cO. 0) OO TO IB 
CS( I ) =CSI 

WLOCxRGASl *CS( I )+RGAS2*( 1 •0-CS( 1 ) ) 

CVM I I ) =CV-1 *CS ( I ) +C V2* ( 1 • 0-CS ( 1 ) ) 

GAML 1(1) =HLOC/CVM( 1 ) 

GAM 1 xGAMLl ( I ) 

I a COnT I NUE 

C.NERC I ) =P< i > / ( 3A.''11*RH0( I ) ) + 0»5 *OVe.L**2 
BOO COnT I NUE 

C ** RETUWN TO REGULAR OPERATION ** 

iKTNxl 


Fig.8-3c (Concluded) 
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** INTtG Ii'iPUT USING USEKIP ♦* 


3 TWO dimensional SCRAMJET SIMULATION *TWO GAS* *SMA«P CORNtW* 


a '2 

.1 000 

lOOO 

200 

0 

-2 

i 

0 

1 0 

1 

940 43 

16 

16 

1 

1 

0 

1 




.2000E-5 

1.0 


1000 







1 .OOOOE-J 

0.0 


1.27 


1 .4 


23.0 

28.0 

- 1 . 0 

.5 
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**** wPoATE i''IuD to 0;»£ 'obERlP* bwAKuOTlNt ♦■»** 

* THPEE Ui.'-'.rivSI 0.-4AL SCRAMJET CAbE * 

\ 

\ 

* I DENT U5E.HIP 
*0 PAGES. 142 

*♦ special. -UStKIP- FOR 3-0 SCrtAiViJtT CASE ** 

OlivIb'NSlON PA.MoC 10) »0AT<5» 10) 
lOi'i Oi 'll 5* 144 ) 

KiiAL) ( b.b'9'0 ) XC.yC 
KEADCbibBO) Jv'lUM 
w)0 ejfeO 

RcAu>( 5*590 ) PANo ( J ) * ( OAT ( SO* -< > » pw= 1*5) 

560 continue 
boO F^RMAT<I5) 
b';^3 T00 .v.aT<6E10,A) 

1 = 1 
I NC = 1 
N-T0T=494 

00 400 tC= I , 1 5 
ITAW= 1 

IF(<.E^.l) ITA,N = 0 
DO 295 i>1=1«nTOT 

1 1 =V./494 
I I =ivi- I I *494 
IF( I I .£0.0 » I i =494 
IIMAX=191 
w/J= I 

IF( I I .LT. I I .-tAX ) GO TO 600 
111=11/19 
111=11-111*19 
1F( I I I .£0.0) 111=19 

I I I VX= 1 O 

1F( I I I .LT. I I li'-iX ) GO TU 730 
Awn= (Y(I)-YC)''(X(I)-XC) 

ANG= ATANC Awg) '*57.295 
JJ= 1 

1F< ANG.L£.PA;'|U( 1 ) ) GO TO 600 
JU=- JNUM 

IF ( ANG.LT.Pa.NG I Ji-lUM) ) GO TO 6bO 

eco continue 

WI =DAT ( 1 * JJ > 
vJl =DAT ( 2* JJ ) 

Vl=OAT(3«JJ) 

Wl=DAT(4.JJ) 

PI=DAT(5« JJ) 

GO TO 750 
650 continue 

Du 700 U=2.jNJM 
uJ= j 

1F( ANG.GT.PaNGI J) ) GO TO 700 
. J1 =UJ- 1 

RRR= ( ANG-PAnG( J1 ) )/(PANG< JJ)-PANG( J1 ) ) 

RI =DAT( 1 • Jl )+XHR*(UAT ( 1 * JJ)-DAT( 1 . J1 ) ) 
ul =DAT< 2* J1 ) +KKR*( DAT ( 2* JJ >-DAT ( 2* Jl ) ) 

Vl=DAT(3-*Jl ) + KRR*(DAT 13 . JJ)-DA7 ( 3 . Jl ) ) 

Fig.8-4c - Three-Dimensional Scramjet Case (USERIP) 
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w1=DAT(<m Jl )+f<rtR*<OAT (4 • )“DAT< 4 « ) » 

P1=0AT<5* J1 )+KRK*(DAT(5« J'J>-DAT(5* Jl » ) 
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7c 0 CONTINUC 
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29b continue 

400 continue 

I WTNx 1 
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Fig. 8-4c (Concluded) 
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Referring to this figure, the nozzle has a sharp 20 deg upward turn at 
X = 0, a sharp 6 deg turn downward at X = X^, and a sharp 6 deg turn outward 
along the dashed line located on the sidewall. The rectangular nozzle exit 
plane is inclined 30 deg from vertical. The body has a sharp 20 deg turn 
upward at X = X^. 

The problem was analyzed in three parts: (1) the nozzle; (2) the external 
flow upstream of the nozzle exit, and (3) the entire flow downstream of the exit 
plane. The nozzle flow field, which is three-dimensional, was computed using 
the GIM code. It was computed using 2160 nodes and required 570 iterations 
to relax to the steady solution. The external flow was computed using a 
Prandtl-Meyer expansion at the 20 deg sharp turn. The nozzle exit plane 
conditions and the freestream expansion were used as input boundary condi- 
tions for the downstream shear layer region. Details are given in Ref. 1. 

The shear layer region is analyzed using 7904 nodes in the GIM code. 
Figure 8 -4b shows the input cards for DYNMAT and GEOM modules for the 
shear layer region. Again the USERIP option of input to the INTEG module 
was used. A FORTRAN listing of this code is given in Fig.8-4c for the three- 
dimensional Scramjet case. The input to INTEG corresponding to this code is 
shown in Fig. 8-4d for the DYNDIM and INTEG modules. 

8.5 THREE-DIMENSIONAL DUCT EXAMPLE CASE 

Next we consider three-dimensional inviscid flow in a square duct as 
shown in Fig. 8 -5a. The first 10 X^ consists of an expansion region with the 
final 10 X being a constant area section. The three-dimensional GEOM and 
DYNMAT deck setups are shown in Fig.8-5b. An example GIMPLT deck is 
shown in Fig. 8-5c. This case is shown as an example of a three-dimensional 
plot setup only; thus the INTEG module is not needed for this illustration. 
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Fig. 8-5b - Three-Dimensional Duct Case (DYNMAT, GEOM) 
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